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William Frederick Durand 
1859-1958 


Ll ribute 


By Theodore von Karman 


| DOUBT that I am the right person to appreciate the 
whole scope of achievement of William Frederick 
Durand and his influence on American engineering 
during the long period of his activities. As a matter 
of fact, I first met him in 1926—two years after his 
official retirement. The fact, however, that I had a 
close connection with him during the last 30 years 
which for anybody else would have been his ‘‘old age’’ 
perhaps enables me to extrapolate as to what must 
have been the characteristics of the earlier periods 0! 
his memorable career. 

Durand was born on March 5, 1859; he received his 
higher education at the U.S. Naval Academy from 
which he was graduated with honors in ISSO. He 
served 7 years in the Navy Engineer Corps. Between 
IS87 and 1924 he was, consecutively, professor at the 
Agricultural and Mechanical College of Michigan 
(4 years), professor of Marine Engineering at Cornell 
University (13 years), and finally professor of Mechan 
ical Engineering at Stanford University. It appears 
that during his period of teaching at Cornell, his main 
interest was concentrated on problems of ship resist 
ance and ship propulsion. With the advent of me 
chanical flight, he shifted his attention and research to 
aerodynamics and flight mechanics. Maybe his great 
hobby, the systematic investigation of aircraft pro 
pellers, derived from his earlier studies of ship propellers 

As had happened to many of the elder members ot 
our engineering community, the First World War and 
the essential role that aviation played in that war had 
a decisive influence on Dr. Durand’s immediate 
interests and activities. 

In 1914 a move was made toward the establishment 
by law of a national aeronautical laboratory. Dr. 
Durand was called to Washington to join in this 
undertaking. When the National Advisory Committee 
for Aeronautics was established in 1915, he became its 
Chairman, through the whole war period. Late in 
1917, he performed a special service as Scientific 
Attaché, at the American Embassy in Paris, in order 
to improve and enlarge the applications of science to 
the military problems of World War I. 

As I write these lines—mostly based on information 
I received recently—the memory of my first visit to 
Stanford returns to mind. Durand was the moving 
spirit of a devoted group of research men and teachers. 





Also, practically every problem, facility, 
was originated by his active spirit which combined 
deep insight with intellectual curiosity. “I am just 
looking for a younger successor,’ he told me. “If 
you decide to come to this country, you can have my 
chair.” I was highly honored, but I was not yet so 
far in my plans, and later Robert A. Millikan’s power 
of persuasion, my own inclinations, and my late 
sister’s liking of Southern California put me in charge 
of the Guggenheim Graduate School at Caltech. 
However, I very soon came into close contact with 
Dr. Durand, through the project of the ‘‘Aerodynamic 
Theory.’ I believe the first suggestion for a hand- 
book that would put down the state of the art as far 
as the scientific fundamentals of aerodynamics were 
between 


every 


concerned originated from a _ conversation 
Daniel Guggenheim and myself. However, Harry F. 
Guggenheim was the person who addressed himself 
with the plan to William Frederick Durand and per- 
suaded him to take over the task of editor. Durand 
conceived the idea having an authorship inter- 
nationally dispersed, over the United States and almost 
all European countries. Thus, he also asked me to 
help him find competent authors who would also be 
willing to cooperate in the undertaking. 

The understanding between Dr. Durand 
Guggenheim Foundation was that the latter would 
cover all honoraria, editorial, and authors’ fees. Never- 
theless— and I consider this as a sign of the change in the 
times between 1930 and today, as far as general interest 
the leading sci- 


of 


and the 


in aviation is concerned—none of 
entific publishers in the United States was willing to 
take the risk of publishing the series of volumes with- 
out a cash guarantee. Even the list of outstanding 
authors, whose names were known over the whole 
world, did not constitute a sufficient security for 
financial success, and Dr. Durand was already willing 
to raise by his own means a sizable sum for the re- 
quired cash guarantee. Then I told him: ‘Please 
give me a few days, so that I may consult my old 
friend Springer in Germany.” A telegram to Springer 
provoked an enthusiastic This explains 
the fact that the six volumes of the ‘‘Aerodynamic 
Theory,” edited by W. F. Durand, and supported by 
the Guggenheim Foundation for the promotion of 
aeronautics, were published in Germany, by Ferdinand 


response. 


Springer, at Leipzig. 

I had the good luck of meeting Dr. Durand on 
many committees and consulting boards. He was 
called upon to preside over the Navy’s special com- 
mittee appointed to study and to make recommenda- 
tions concerning the future of lighter-than-air air- 
craft. I was also a member of this Committee, prob- 
ably on account of my previous experience as con- 


sultant to the Zeppelin Works in Germany. Dr. 
Durand and I were also associated with the U.S. 


Bureau of Reclamation, where he participated in the 
planning of the great dams and hydraulic projects 
in the broadest sense, while I worked with the late 
Robert T. Knapp and with A. Hollander on the more 
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During th 


restricted aspects of pump development. 
Second World War, when Dr. Durand, at the age of §9 
yas called upon to become Chairman of the specia| 
NACA committee on jet propulsion, our interests me 
in this most recent field of aerial propulsion. 

I believe that the main characteristics of Dr. Dy 


rand’s influence on American engineering can 
expressed in a few words. American engineering 
of the Nineteenth Century suffered—ike engineering 


of all countries—from an overestimation of the 59 
called practical engineering and an underestimation 
of the possibilities of theoretical insight and theoretica| 
computation. The main gap existed between hydrau. 
lics and rational fluid 
process between experimental evidence and_ theory 
was started in England by men like Osborne Reynolds 
and Lord Rayleigh; in Germany, to some extent by 
Helmholtz, and certainly by Felix Klein and Ludwig 
Prandtl; whereas in France, where the purely mathe- 
matical theory of fluid mechanics was especially well 
developed, the gap survived until a late date—possibly 
because of the traditional influence of the Ecole Poly- 
technique. In the United States, W. F. Durand 
carried the torch for the approach between theoretical 
Owing to the great repu- 


mechanics. The conciliation 


and practical engineering. 
tation he enjoyed among the followers of traditional 


empirical engineering, he was able to educate a genera- 

° ‘ ‘ . 7 

tion to the modern views on theoretical treatment of 
t 


engineering problems. I do not believe it is an exag- 
geration if I ascribe the most important role in this 
change of views and methods to Durand’s activities 
In a way I always felt that my own work—after | 
came to this country—was greatly facilitated by his 
influence on the American scene of the engineering 
science, and I am sure that my friend S. Timoshenko, 





who fulfilled an important analogous role in the field | 


of elasticity and structures, will confirm my judgment. 

I met Dr. Durand for the last time in 1955, when I 
gave a lecture on ‘“‘Solved and Unsolved Problems of 
Modern Aerodynamics” on the occasion of the cen- 
tennial celebration at Brooklyn Polytechnic Institute. 
He expressed the desire to attend my lecture. Since 
the death of his beloved wife, in 1952, he lived in his 
son’s house in Brooklyn, and Dr. Hoff, then Head of 
the Aeronautics Department of the Institute, sent a 
car for him. After the lecture, Dr. Durand told me: 
“You know, Karman, in the last one or two years 
I really began to feel that I am getting old’’; he was 96 
that year. He suffered some deterioration of his 
eyesight but was otherwise physically and mentally 
alert. 

Unfortunately a stroke came about a year later, and 
our hope that he would be physically present at a 
Durand centenary celebration was frustrated by his 
passing away on August 10, 1958. 

At least two generations of 
and many engineers and scientists from all countries, 
who had the good fortune to meet and know him, 


American engineers, 


will long remember him and cherish many happy | 


recollections. 
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A Theoretical Investigation 
Between Shock Waves and 


otf the Interaction 
Boundary Lavers 


M. HONDA* 


Tohoku Unwersity, Sendai, 


SUMMARY 


This paper treats of a self-inducing type of flow which appears 
in the interaction process between shock waves and laminar and 
turbulent boundary layers. In the case of laminar boundary 
layer, the development of the layer is described on the usual 


assumptions. In the case of turbulent boundary layer, the layer 


is divided into two regions 
to the wall and the other the outer inviscid layer. 
growth of the inner viscous layer following Karman’s dynamic 


-one the inner viscous layer adjacent 
Then the 


similarity law is determined by the flow conditions at the inner 
edge of the outer inviscid layer. The obtained results agree 


fairly well with those of the experiments 


SYMBOLS 


coordinates paraliel and normal to surface 


x, ¥ = 

u, = components of velocity in x and y directions 
6 = flow direction 

6 = thickness of boundary layer 


6* = displacement thickness defined by 


é 
f [1 — (pu/pim,)|dy 
0 , 


b = momentum thickness defined by 

i) 

(pu/pitts:) [1 — (u/m,)|dy 

o* = energy thickness defined by 

Cr) 

f, (pu/pits) [1 — (u?/u,?)|dy 

P, p pressure, density 
¥ = ratio of specific heats 
‘fs = absolute temperature 
a = velocity of sound 
M = Mach Number 
m = [(y — 1)/2}] 1? 
#,v = viscosity, kinematic viscosity 
T = shear stress 
X, Y = coordinates in the transformed flow, Eq. (1.4 
H = }[(6*/3) — m]/(1 + m)}, Eqs. (1.7) and (2.13 
G = 9*/9 
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Japan 
lag = form parameter for velocity profile in turbulent 
boundary-layer flow, Eq. (2.14) 
d; = turbulent boundary-layer thickness defined by 
, (u u,) [1 (u/u,)|dy, Eq. (2.10 
U; = friction velocity, (ry/ pu yi/e 
U, U, = u, u; in the undisturbed flow 
f = universal function for velocity profile in turbulent 
boundary-layer flow, Eq. (3.8) 
y¥ = stream function 
¥r = U;zy/r 
y* = (U, re) fy 0 pwidy, Eq. (3.7) 
Subscripts 
0 = conditions in the free stream 
1 = conditions at the outer edge of the boundary layer 
w = conditions at the wall 
i = quantities transformed by Eq. (1.4) and by Eq. (2.11 
s = conditions at the point of the laminar boundary-layer 
separation 
d = conditions in the dead air region 
eé = conditions at the inner edge of the outer inviscid layer 


INTRODUCTION 


I MANY AERONAUTICAL PROBLEMS—such as the flow 


N 

I past wings and bodies, the flow through engine 
intakes, rocket nozzles, and wind-tunnel diffusers 

shock waves interact with the boundary layer on the 
surface. In such circumstances, the pressure distribu- 
tion on the surface has no discontinuous jump, but a 
continuous rise. Particularly, in case where the ex- 
ternal flow is of the simple wave type, the pressure dis- 
tribution upstream of the foot of shock wave depends 
mainly on the free-stream Mach Number and the 
boundary-layer Reynolds Number and hardly at all 
on the agency which provokes shock wave. This fact 
explains the existence of a flow to be determined by an 
equilibrium between two processes: the pressure rise 
is induced by the external flow deflection brought about 
by the thickening of the boundary layer, and the thick- 


ening is governed by the pressure rise. The present 
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paper treats of such a self-inducing type of flow in cases 
of laminar and turbulent boundary layers. 

In Section (1), the development of the laminar 
boundary layer under consideration is described by 
using the momentum- and energy-integral equations 
on the usual boundary-layer assumptions. The ob- 
tained results are compared with the Gadd! and the 
Gadd, Holder and Regan experiments.” In Section 
(2), the development of the turbulent boundary layer 
with adverse pressure gradients in incompressible and 
compressible fluid flow is treated semiempirically on the 
usual boundary-layer Then the ob- 
tained calculating formulas are applied to the approach 
to the subject, and it is pointed out that a theoretical 
contradiction lies in such a course of approach. It 
seems to the author that this shows the inadequacy of 
applying the development of the turbulent boundary 
layer as deduced from the one in incompressible fluid 
flow, as it is, to the flow with abrupt pressure rise as 
brought about by the interaction with shock waves. 

In Section (3), it is assumed that the turbulent 
boundary layer under consideration can be divided into 


assumptions. 


two regions—one the inner viscous layer adjacent to the 
wall and the other the outer inviscid layer as in Light- 
hill’s theory.’ The growth of the inner viscous layer 
following Karman’s dynamic similarity law‘ is then de- 
termined by the flow conditions at the inner edge of the 
outer inviscid layer. The pressure distributions at the 
wall obtained thus are in fair agreement with the Gadd! 
and the Fage and Sargent’ experiments up to the point 
of the boundary-layer separation. 


SHOCK-WAVE - LAMINAR BOUNDARY-LAYER 
INTERACTION 


(1) 


Basic Equations 
Consider first the steady two-dimensional boundary- 


layer flow on the usual boundary-layer assumptions. 
Integrating the continuity equation with respect to y 


6,/(1 + mo) = 3(dH/dX) + )H + [m,/(1 + my) |} (dd,/dX) + 


(cay — 1)/(v¥ —- )DJA+ 24+ [(¥ + I/(7 - 


(apM,/v) (dd ?/dX) + 2 (2 + H) X 


(aod 7/9) (dA, /dX) = PU) (1.9) 


(ay MM, Yo) (d dX) (Gd ;*) at 


6(apG?d ;?/ v9) (dA, /dX) QUT) (1.10) 


where 
P = 2(9,/m, 1) (Ou; OF), 
QO = 4(Gd;/m, :”) { (0u,/OY)*dY, G = #,*/9, 
si (1.11) 


Moreover, eliminating dd3;/dx from Eqs. (1.9) and 
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from y 0 to y 6, we have the following mass-fly Mi 
equation 
= . 1) Mie 
d6* 1 d(pyt1) "pu ; 
6; — dy ] ] 

dx pity dx Jo ply where 
where 6; is assumed to be sufficiently small so that tay Pia | 
6, = 0, and cos 6; 1. The momentum-integral equa 
tion is 

vi 
(1 / pitty?) (d/dx) (pity?) + 
(§* ( - T 2 9) ; - 
6* /u1) (du,/dx) » prlty 1.2 Consi 
Provided that the zero heat transfer at the wall and 
7 ° ° w0° 
the Prandtl Number of unity are assumed, the energy 
integral equation can be written’ 
3 ] pe 3.9 * Q * 2 2 ] hold i} 
(1 / pytty®) (d dx) (pytty?d3*) + 9 (* / 017) (day? da 
2 | (7 pitts”) (OO) (u/m)dy (13 aC 
J0 
The three equations above hold in the cases of both 
laminar and turbulent boundary layers. it 
. : ind ti 
Development of the Laminar Boundary Layer 
For the laminar boundary layer, on the assumption 
that the viscosity is proportional to the absolute tem hen 
le a 


fluid can be transformed into one in incompressible fluid 


perature, a laminar boundary-layer flow in compressible 
by the well-known Stewartson transformation® 
| 





, , j 
dX (Pid Pod dx, dy (pid, poly) X ( @ 
(p/pi)dy (1.4 
M1, ay Vi, Ui /M, § = U/ ty 
é* (potty prdi) }(1 + my) 6;* + my (15 
I = (poly pidi)d;, d3* ( polo pid, )d;™ (1.6 dH 
Assuming that the velocity profile in the transformed 
layer can be expressed as a member of a one-parameter 
family of curves, we introduce the form parameter 
. 7 Howe 
H 6;*/0 (Li 0(10 
By virtue of Eqs. (1.4), (1.5), and (1.6), Eqs. (1.1), | PUTPS 
(1.2), and (1.3) can be rewritten Eqs. 
[m2,/(1 + m)| X 
1)] [aay/(1 + my)| + [(MQ? — 1)/m(1 + m)| X 
} (u;/1, ;) (dV a) (3;/M,) (dM,/dX) (1S D 
J 0 
Lami 
(1.10), we have 
Wi 
(1/G*) (dG? dH) [(ao\h/ )d ?(dH dX) | of tk 
[(Q/G?) — P| + 24 — 1) X to an 
[(aod 2? vo) (dM, dX)|] (1-12) ) mutt 
The relations between //, P, Q, G, etc., are deter- | wall 
mined conveniently by adopting the Pohlhausen’s oer 
fourth degree polynomial’ for the velocity profile in on 
the present paper. Table 1 shows their numerical fore, 
values. stam 
f the v 


Inserting Eq. (1.9) into Eq. (1.8), 


SS-fly) 


ll and 
eT LY 








«| 
(13 


” 


both 


potion 
tem 

sible 
fluid 





1S 


INTERACTION BETWEEN 


Mads v0) [1/(1 + m0) 
} H+ [my,/(1 + m,)]f (P 2) + 
Md; vo) (dT dX) — Flag? vo) (dM, dX) 1.13) 
where 
F= H(2 1) — [2y/(4 1)| [vy m,) ZT — 
y+1)(y — 1)] [nn/ 1 + mm) - 
V/ 1) (1 + m,)? (u; 1, 3) (€Y 9;) 1.14) 
J0 
Considering that the relations 
3)" Po(voX /Uo), (Uo8o/v0) [00/(1 + mo) 
} Hy + [mo/(1 + mo) ]f (Po/2) (1.15) 
hold in the shock-free flow state, we rewrite 
8, /(1 + mo) = (v0/ Udo) } Ho + [mo/(1 + mo) ]f x | 
(P,/2) + [0,/(1 Mo) (1.16) 
o,? Po(veX / Uo) (8/80)? 
ind take the variables as 
d= dX/X, 19 In (14o/M,), ¢ In (3; BJ) (1.17) 
Then Eqs. (1.9), (1.12), and (1.13) are 
dt dé) (2 T IT) (dn dé) 
—(1/2) [1 — (P/Py)e’~* 1.18) 
d In G dil dn 
- — (71-1) 
dH dé dé 
1 /O ‘ 
-——<—- re 1.19) 
2Po \G? 
dH /dé) + F(dn/dé) (Updo/ volo) X 
[0,/(1 + mo) lef + (1/2) GH 4 
[m9 /(1 + mo) je '-j}H + 
[m,/(1 + m)]} (P/Poe’~*) (1.20) 
However, from the fact that &, ¢ O(10~!) and 


0(10~*) in the actual phenomena, let us put, for the 
purpose of approximation, é, e”, ¢' 1, my my in 
Eqs. (1.18), (1.19), and (1.20), and 


SHOCK 


W 








AVES AND BOUNDARY LAYERS 669 
PABLeE | 
| 
- din & Q/iG?*)-P \ ° sy | 
" rey rw aPUF 2 P 
2.6 0380 0014 4.973 5 .55¢ 
ee ~0295 OC 4205 5.666 
2.0 022 0725 4500 5 
2.9 0174 1 2554 
3 014¢ 1226 2263 47 
3.1 0094 141 -172yu 421 
4.2 64 15604 2 4 5.207 
3.43 .0036 1606 0734 5.2606 
3.4 0018 1786 037k 5.254 
3.5 000! L664 000 5 .25C 
H 2.597, Po = 0.4393; H, = 3.500, G 1.5315, (O/ 
0.4004 
ISENTROPIC FLOW 
2 
SEPARATED LAYER 
u=O LINE 
ee : 
es 5, DEAD AIR REGION 
s WALL 
Fic. 1. A model for the flow downstream of the point of laminar 
boundary-layer separation 
Updo/voP V Ro/V Po, R Uoxo/v (1.22) 
where x» is the distance from the leading edge of the 


plate to a point which represents the phenomena.t 
Then we have the following calculating equations for 
bond > | 


0,/(1 + mo) [V Mo? — 1/(1 + my)? |n (1.21) 
obtaining the pressure distribution at the wall 
dé/dH = (1D) j\1 — [F/(H — 1)| (diInG di); 
dn/dH = —(|1/(H — 1)] (dInG/dH) + }[(Q/G*) — P| 2P,)U7 — 1); (dé di) | 198 
.23) 
D = [VU Mo? — 1/(1 + mo)?] (4/ Ro V Po) 9 + (1/2) Go + [mo/(1 + m) |} - \ 
(P Po) |} H + [mo/( + mo) ]t — [F/ Po( 1)} (0G?) P}) 
Laminar Separated Layer 
We now consider the flow downstream of the point 
of the laminar boundary-layer separation. In order Consider a flow model as shown in Fig. 1. Assuming 


to analyze such a flow precisely, we need to clarify the 
mutual relations between the layer separated from the 
wall and the dead air region adjacent to the wall. How- 
ever, it seems next to impossible to solve the Navier- 
Stokes equations in detail at the present stage. There- 
lore, in approaching the subject under such circum- 
stances, it is necessary to set up some assumptions on 
the way. 


that the usual boundary-layer assumptions still hold 
in the separated layer, we have the following mass, 
momentum, and energy-flux equations for the sepa 


rated layer: 


of the 


nomena under consideration to be small compared with the dis- 


+ This approximation expects a priori the scale phe- 


1 


tance from the leading edge. 
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pit, ((d6 dx) + (dba dx) — 6; | + pda = | 
ba+s . Ls} 
(d ax) { pudy (1.24) 
bd Pp 
Po 
(1/ pyr?) (d/dx) (pytti?d) + (6*/u1) (du, /dx) 
Pata pilty (1.25) 12 | . 
(1 pitt”) (d dx) ( pyut;>3*) + m(o* U7) (du;" dx) = Fé 
bats 4 
(pata pitti) + 2 (r/ pitts”) (0/Oy) (u/ui)dy (1.26) 7% ° GADD S EXPERIMENT 
ba / - THEORY 
J. 
Eq. (1.24) can be rewritten as ft 
aa lof 
6, — (déy dx) = (dé* dx) — (1/ pitts) [d( piu) dx] X -0.2 “61 0.0 0.1 02,703 
dbd+sd ° 
= x ol 
f (pu! pitti)dy + (pava pitti) (1.27) 
F) Fic. 2. Pressure distribution at the wall when JJ, = 3 
and Ry ~ 4.1 K 10°. © = interaction with the impinging shock 
wave generated by a wedge in the free stream; XK = interaction 


By using the Stewartson transformation, we have 


with the shock wave generated by a spoiler on the wall 


[1 (1 + mo)] (0, — (dba dx)| = 3(dH dX) + ;H + [m,/(1 + m)]} (dd;,/dX) + [m, (1 + m)] X 
By -— 1) (vy — DIA +2 + My + D/(y — ‘DY [me + om) | + OC? — 1)/m(1 + m)) X 


f (uw; (m3, 1.) (€V/3,)} (8; My) (dM, /dX) + [1/0 + m)] (vam); (1.28 


(dd; dX) + UT + 2) (8; Mh) (dM /dX) = 


(Va 4) (1.29) 
(dd ,*/dX) + 3(8,;* M,) (dM,/dX) = 
2(v 14, »f (Ou;/OYV)*dY + (vg/u1); (1.30) 
where Va, i = (pa/ pi) (Ao/a1)* va 


We now set up two assumptions: one that the inner 
edge of the separated layer is linear, and the other that 
the velocity profile in the layer is same as the one at 


the point of the boundary-layer separation. That is, 


(F — {(H, — 1) [1 — (1/G),]}) (dn/dt) + (X/X, 


v;) (dl dé) 


assume a jet-like flow, then 
[1/(1 + mo)] (déa/dx) = [8;(dI7/dX) | (1.31 


in Eq. (1.28), where subscript —s denotes the value on 
the upstream side of the separation point, and Eq. 
(1.30) is 


(dd; dX) + 3(8;/M,) (dl, /dX) = 
(Vo aoM,d;) (OQ 2G?),. +- (1 G,) (Va ,) (1.32 


Under the same procedures as in the case of the attached 
layer, we arrive at the following equation: 


,+ (U1, + 1)/2Po] Xx 


1(Q G?), (1 = (1 G),}} e" = (Udo voPo) [O; fq. + my) | e . + (1 2) ) Ho + [m1 (1 + mo) |te . (1.33) 


Moreover, on the same approximations as used in obtaining Eq. (1.25), Eq. (1.33) is simplified as 


(2 — 1)/{1 — (1/G),]} — F) {(dn/dt), — (dn/dé)} 
Solving this equation, we have 


n — ns = [(Q/G*),/2Po(, 


«= (WM? — 1/(1 + m)?] (WRo/V Po) CG (2 — 1)/[1 — (1/6), |} -— FDS 


with respect to the pressure distribution at the wall 
down the point of the laminar boundary-layer sepa- 
ration. 


Comparison With Experiments 


We now attempt to compare the results obtained in 
the theory with the experimental data. The experi- 
mental measurements of the pressure distribution at 
the wall by Gadd! under the interaction with the shock 
waves (one the impinging wave generated by a wedge 
in the free stream and the other the wave generated by 
a spoiler on the wall, when J) ~ 3.0), are shown in 


[V My? — 1/(1 + m)?] (V Ry V P,) (n — n;) (1.34) 


1)} (1/«) (A — e7*47*) } 


Fig. 2 in comparison with the theoretical curve. The 
origin of the x coordinate is positioned at the point 
where the separation occurs according to the theory, 
and the Reynolds Number depending on the distance 
from the leading edge to that point is 4.1 K 10°. The 
measured pressures increase sharply at some distance 
downstream from the separation point because the 
transition to turbulent flow takes place there. How- 
ever, Fig. 2 predicts that the pressure at the transition 
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pressure curve—and the terminal pressure rise of the 
laminar separated layer according to the theory are 
the quantities comparable with each other, because 
the pressure at the transition point may scarcely de- 
pend on that position from the characteristic of the 
saturation-type curve. 

Let us now compare the terminal pressure rise in the 
laminar separated layer according to the theory with 
the transition pressure measured by Gadd, Holder, and 
Regan’ in the cases of the interaction with an impinging 
shock wave and a shock wave generated by a wedge on 
Fig. 3 shows the comparisons 
The results 


the plate, respectively. 
between them when J/) ~ 2.0 and 3.0. 
obtained in the theory agree fairly well with the experi- 
mental data for the simplified analysis in it. 


(2) DEVELOPMENT OF THE TURBULENT BOUNDARY 
LAYER WITH ADVERSE PRESSURE GRADIENTS 


As mentioned in Section (1), the interaction process 
between the shock wave and the laminar boundary 
layer can be described approximately by combining 
the development of the laminar boundary layer on the 
usual boundary-layer assumptions with the relation 
between the magnitude of velocity and the deflection 
angle of velocity from the free-stream direction at the 
outer edge of the boundary layer. Let us attempt to 
approach the case of turbulent boundary layer along 
similar lines, putting aside, at present, the question of 
adequacy of this approach. In such a course, the first 
question is to describe the development of the turbulent 
boundary layer by using some equations. In the pres- 
ent stage where there is no experimental data as there 
is for the development of the compressible turbulent 
boundary layer with adverse pressure gradients (in 
particular, the change of the velocity profile in the 
layer), we are compelled to use the analogy of the de- 
velopment of the layer with adverse pressure gradients 
in incompressible fluid flow and the one with zero 
pressure gradient in compressible fluid flow. 


Incompressible Turbulent Boundary Layer 


A number of calculating equations for this subject 
have been put forward already, and the main difference 
between them is in the expression for the nondimen- 
sional work done by the shear force, 


E= [ (r/ pur”) (0/ Oy) (u/m)dy (2.1) 
J 0 


which appears in the energy-integral equation [Eq. 
(1.3)|. The expressions about this term, or the term 
including this, presented by Gruschwitz,'’ Doenhoff 
and Tetervin,'! Garner,'* Truckenbrodt,’ Rubert and 
Persh,!* Schuh,'4 and Tani™ are considerably different 
from each other, and it is questionable as to which 
should be adopted. Therefore, in the present paper, we 
derive the expression for the term 


IT = [1/(4 — 1)] [2(E/G) — (tw/pm?)] (2.2) 


from the experimental data by Gruschwitz," Doen- 
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hoff and Tetervin,'! and Schubauer and Klebanoff.'* 
In the following equations for the incompressible 
turbulent boundary layer 


(dd dx) + (+ 2) (8/m) (du,/dx) = 7, puy-t a 
(dd* (dx) + 3(8*/u,) (du,/dx) 2E f Pap 


G = 3* J could be approximately a function of // alone, 
and there holds Wieghardt’s semiempirical formula” 


G = 1.269H (HW — 0.379) (2.4) 


Eliminating dd dx from Eqs. (2.3) in virtue of Eq. (2.4), 
and integrating it with respect to +, we have 


| I(dx/d) In [(t0/21) (Lo/L) (2.5) 


where x is an arbitrary point and 
L = G{(H — 1) (H — 0.379) |" 


Then, in the cases when the velocity, //, and J dis 
tributions are presented, we can calculate numerically 
the term in the right-hand side of Eq. (2.5) and obtain 
the numerical values of / by differentiating that term 
with respect to x. It should be pointed out that the 
data used here are restricted to those in the domain of 
the so-called fully developed turbulent flow —i.e., the 
domain dH] dx > 0. 


Tillmann’s experimental formula’* for the wall shear 


In consideration of Ludwieg and 


stress, 
A(H)/R;, Rao ue v (2.6) 
A(/T) 


Ty ply” 


where ¢ = 0.268, 0.123 & 10-0-68# 


we now arrange R,’/J obtained above against a param 
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Fic. 5. Wall shear stress against the Reynolds Number depend- 
ing on the thickness J; defined by Eq. (2.10). 
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eter Ry’ °(0 uw) (du, dx) and then have an approximate 
relation 


RJT = 0.0018 — 1.08R,° (8 /u;) (du, dx), 
(1.3 << H < 1.8) 


(234) 


Consider briefly the foundation 
According to Fediaevsky, 


as shown in Fig. 4. 
on which Eq. (2.7) stands. 
the distribution of the shear stress near the wall is 


T Tw + (dpi/dx)y + O(y’) 


If we adopt the form of the velocity profile as 


Ou/Oyx u,/¥, U, W Fu/2 


a priori that the flow near the wall contributes 


expecting 
i.e., to the 


largely to the work done by the shear stress 
creation of turbulence——a functional relation about the 


work done by shear stress could be assumed as 


f£ = function of 
[7./ pl", (6/ pty") (dp, /dx)(u,/m), H| (2.8) 
It is clear that Eq. (2.7) is a simplified expression of 


those deduced from Eq. (2.8). 


Compressible Turbulent Boundary Layer 

What kind of form is to be assumed of the velocity 
profile in the compressible boundary layer? It is 
questionable. It has been shown that a number of 
measured velocity profiles in the compressible turbulent 
boundary layer with zero pressure gradient could be 


roughly represented by the power law 
U/Uo = (y/6) 


and Karman’s dynamic similarity’ holds closely to the 


wall, 1.e., 
U/U, = f(U,y/v«) 
Then, in the present paper, we assume the following 
form: 
u/u, = (y/6)'" y>vd) 
: , oe (2.9) 
uU/u, = J(ULy/Vo) VY > VB; 
6 
where o, = [ (u/u1) [1 — (u/m) |dy (2.10) 
/7 0 


just as in incompressible turbulent boundary layer.*° 
Now, introducing the transformation 


I = (podo/pia))d;, 3* (poo / pid) 0 ;* = (2.11) 
in Eqs. (1.2) and (1.3), we have 
(podo/ puti) [(dd,/dx) + (H+ 2) X 
(3;/M,) (dM,/dx)| = ty vst} 
(poo / Pots) [(dd;* /dx) + 3(3;*/M)) (dM, /dx) | =) (2.12) 
2 | (7) pwtti?) (O/Oy) (u/m1)dy 
JV 
where H = [(6*/3) — m]/(1 + m) (2.13) 


Tables 2 and 3 show the relation between the equation 
parameter H7 and the form parameter for velocity pro- 
file 7* as 
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Eg 1+ (2 n) 2.14) 


and the one between the momentum thickness J and 

the thickness 3; defined by Eq. (2.10), respectively. 
From Eqs. (2.9) and (2.10), we could deduce that 

the following functional relation about the nondimen 


sional shear stress at the wall holds: 


Tw / Ptr" function of (m0, v 


Particularly in the flow with zero pressure gradient 


Tw Poo” must be a function of Ly, v, alone because 
I1T* would be referred to ly, v,, in this case. 


this inference, the experimental data by Coles*' are 


To justify 


rearranged in Fig. 5 on the assumptions of w= 7° (w 


0.76), recovery factor O.SS and n 7. Thus, we 
have 
Tw/ Potts" A (F1*) (ud; v 2.16) 
corresponding to Eq. (2.6) and 
‘ I ta jv Ere 
(10;3;/V.) [1/(* — 1) ) I= G(H*)| X 
{ 
(7 Plt) (O Oy) (u m)dy Tw Plt”) ¢ 
J0 
0.0018 1.OS( 14,2; v,,) . (py Pw) (9; /,) (diy dx) 
(2.17) 
in accordance with Eqs. (2.7) and (2.8). 
Consider the relation between G and // or J/*. We 


estimate the values of G approximately by substituting 
the power law for velocity profile in G and arrive at 
the result of the compressibility effect on the relation 


between G and // small enough up to 1/; for 5, as 





1.8 
7 
M -5 
M,= 3 
M,= 0 
164+ . nm - 
1.2 14 1.6 1.8 
n 
Fic. 6. Compressibility effect on the relation between G and H, 
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shown in Fig. 6. Therefore, in the present paper, we 
assuine the same relation as the one in incompressible 
fluid flow —i.e., 

G = 1.269H/(H — 0.379) (2.18) 


where // is the equation parameter defined by Eq. 
(2.13). 
Let us derive the equations to describe the growth of 


SCIENCES NOVEMBER, 


1958 


vg and //, respectively, by using the above-mentioned 


equations. Considering that 


U0; Ve = (Ao/Q1)™ (pw! pr) (putt; / pid) (Ao Md 


we carry out the transformation 


y—1)] +29 


dX (dy Q) 


1 
(Pw/ Pi) 


7 


x 


(pit / pyt;)" dx 


into Eqs. (2.12), then have 


(dd ;/dX) + (IT + 2) (8:/M)) (dM, dX) A (a9M,0 ;/%)~ 


and 


1 i 


where B 


on the approximation of identifying G(/7*) and G(/7). Table 4 shows the numerical values of B, where (a; /ao)”’ 


and recovery factor = 1 are taken approximately. 


Interaction With Shock Wave 


We now attempt to apply Eqs. (2.20) and (2.21) to the subject under consideration. 


Eq. (1.1) can be rewritten as 


‘{—1 (7 — 1)] (dInG dH)| 8(dil/dX) + [1 — 1.08B(apAQ8; v9) ~' 7 (3;/.,) (dM,/dX) 


—0.0018 [7* — 1)/(H — 1)] (ao 448; 0) | 


(a;/ao)” (p1/ pu)” . [(p1/ pw) (1 + my)| [(U7* — 1)/(H — 1)] (pid; pi)! 


=) 


(ao/a:)"" (p1/pe)' (pud;/ pd)” [0:/(1 + m)) = 38(dTT/dX) + SIT + [m,/(1 + m)]f (dd; /dX) + 


[m,/(1 + m)]9 [(34 —1)/(y¥ —- I) JH+2 4+ [(y + 1)/(y — 1)] [m./0 + m))] + 


( 


%0 


[((? — 1)/m(1 + mdi (pu/piu,) (dy at (3;/My) (dM, /dX) 


We ignore the terms in the right-hand sides of Eqs. 
(2.20) and (2.21), in consideration of their orders being 
10-* ~ 10~4, and, on the other hand, the scale of the 
phenomena brought about by the interaction with 
shock waves being of comparable order with the bound- 
ary-layer thickness. Moreover, using the approxima- 
tion py / py, 1 + min Eq. (2.22), and taking the vari- 
ables as 


dE = (py; pit)’ (daX/3%), 9 = In (Mo/M), 
¢ = In (8,/80) (2.23) 


we have 


dt/dH = [(1 + m)!*” VM? — 1] (e /n)Jt (2.24) 
J=1-—(1/C) [F/\(7 —1)|)(dinG/dH) f ~~ 
dy/dH = —(1/C) (WT - NI MG oy. 
d¢/dH = (H + 2) (dn/dH) , oe 
where 


C = 1 — 1.08B(ayM,9;/1)~"’? 

F = H(2 + H) — [27/(v — 1)] X 
[my /(1 + m) WT — [(y + 1)/(v — 1)] X (2.26) 
[m,/(1 + m)]? — [(A72 — 1)/(1 + ™)?] X iit 


i} (pu / pitt) (dy 8) 


The question arising in these equations is the sign of / 
in Eq. (2.24). The sign of J is not always plus. In 
fact, as shown in Fig. 7, the domain of J < 0 appears 
as soon as the Mach Number exceeds unity, just as 
the domain called by the word “‘supercritical’’ appears 


7 


Vy) 


2.19) 


(2.20) 


(2.21 


The deflection angle 4, in 


(2.223 

















TABLE 4 
Coefficient B in Eq. (2.21) 

fax} 8 7 6 | 5 | 4 3 | 2.5 
0.0 | 1,000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 
1.0 943 | .949 | .952 | .957 | .968 | .979 | .985 
1.2 922 | 931 938 | .9h6 -959 974 . 982 
1.4 900 .908 920 | 933 -948 -966 977 
1.6 879 | .888 g02 | .919 -937 961 974 
1.8 -855 | .870 886 | 904 -926 953 -970 
2.0 836 | .851 870 | 891 -916 “949 .967 
2.2 -817 | .835 .855 | 880 -908 944 966 
2.4 .800 | -818 82 .870 903 -94.0 -962 
2.6 785 | .806 | .830 | .859 | .89h .936 | .961 
2.8 -772 +793 .819 .851 .887 934 -961 
3.0 -759 | .783 .810 By .882 -932 961 
3.2 748 | 2773 .802 .838 .878 -930 -961 
3.4 -738 | -762 -794 -832 .875 -930 -962 
3.6 727 | 755 .786 .826 .872 -929 -962 
3.8 .718 | .748 | .780 | .821 | .870 | .929 | .964 
4.0 no | .7y2 | .775 | .816 | .868 | .929 | .966 
4.2 -703 | .736 771 -813 | .867 -930 968 
4.4 -697 730 -767 .811 -866 | .931 971 
4.6 -693 | -726 -764 .808 -865 +933 -973 

| 4.8 . 688 -722 -761 -806 864 -934 | -976 

| 5.0 68 -718 -758 80k B64 935 | -979 | 
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domain, 


in the Crocco mixing theory.** In that 
these equations give a solution with no meaning in 
It seems to the author that this contradiction 


reality. 
we apply the 


is due to the following assumption: 
ysual development of the incompressible boundary 
laver where the shear force acts as standing comparison 
with the inertia and pressure forces throughout the 
laver, as it is, to the flow with abrupt pressure rise as 
brought about by the interaction with shock waves 
in other words, to the flow where the inertia and pres- 
sure forces may be superior to the shear force in the 
larger part of the layer. 


(3) SHocK-WAVE - TURBULENT BOUNDARY-LAYER 
INTERACTION 


In taking the superiority of the inertia and pressure 
forces to the shear force into consideration, the course 
of approach founded on the inviscid solution which has 
been developed by Howarth,** Tsien and Finston,** 
Robinson,” and, particularly, Lighthill® is certainly a 
useful guide. However, it seems to the author that the 
question arising in their theories is their treatment of 
the so-called inner viscous layer close to the wall where 
the shear force may act as standing comparison with 
the inertia and pressure forces and, in particular, their 
treatment of the nonlinear phenomena in that layer. 
Therefore, we develop the theory focusing our attention 


on the nonlinear phenomena. 


Outer Inviscid Layer 
Neglecting the viscosity and the heat conduction 
that is, assuming the constant specific entropy along a 


streamline —we have 


(1 — AW?) M*| (0/0s) In p” = 06 On| 


(1/M?) (0/0n) In p’7 = —(00/0s) f me) 


where ds and dn are the line elements tangential and 
normal to the streamline, respectively. As an approxi- 
mate solution, 


In pi? In p,\/7 + (dO; ax) | M? (dy pu) +... 


6 = 6 + [d(In p,'/7) dv] f [(M? — 1)/M?] XK 
(dy/pu) +... 
(3.2) 


be- 
the 


Substituting the following approximate relation 
tween the deflection angle and the pressure at 
outer edge of the layer, 


1)/yl 


0; (VM — 1 2m) }1 — (po/ pr)” ; Ze) 


into Eqs. (3.2), we have the relation between them at 
the inner edge of the layer 


G& = (V M.? — 1/2m)A — [(1 + mo)?/2mo] X 
(1 — AjO/?IB—v/G—-DI Fray, A) (dA/dx) (3.4) 
where A = 1 — (fo/p,)'%~ (3.5) 
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20} 
Re o* 
M R,*10° 
J<O 
St —- 
—_ ron 
ocumaias 14 i6 i8 
H 
Fic. 7. Sign of J in Eq. (2.24 
Fy = [Mo?/(1 + mo)?21 (11 — A)7 OV 1G 41 


‘| ((1 — M2) M2] (dw / pu) + [2 — 1) Mo 


' 
mae ' 
| AM? (dy pu) ¢ 


need to assume 


3.6) 


In order to calculate Eq. (3.6), we 


the form of the velocity profile in the layer. For two 
reasons, one that the region close to the wall including 
the laminar sub-layer plays a predominant role and the 
other that this is the convenient form for calculations, 
we assume the following form” as the velocity profile 
in the undisturbed flow state: 





U/l rs. ¥* (U,/vg) | (p, pr)dy (3.7) 
where / is a universal function written as 
fe pe <1 
J - e4ket - " (3.8) 
2.5in ¥* +55 Y* > 11.79 
Assuming constant stagnation enthalpy all over the 
field, 
p/p = (po'!” po) | 1 + moll — (U? Uv)f IY 4 
1 ¥ * ) 
(1/2)u? + [v/(v — 1)] (p/p) constant f P 
2S} 
LA 
Pe 
v4 
2.0} 
1s} GADD’S EXPERIMENT 
— THEORY 
a ee = 
-10 0.0 1.0 Eo . 3.0 
Ort y i9-* 
~ a 
Fic. 8. Pressure distribution at the wall when My ™ 3.0 and 


Rs ~ 7.5 X 104 
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F, (uw/ poo) (Uo? / U,?) Fe 


where 


Inner Viscous Layer 


Assuming that the usual boundary-layer assumptions 
hold and the compressibility of fluid could be ignored 
in the inner viscous layer, we have 


Pytt(Ou/Ox) + pyv(Ou/Oy) = 


—(dp./dx) + (O7r/Oy), pu constant (3.11) 


(Ow/Ox) + (Ov/Oy) = 0 (3.12) 
Moreover, we assume the velocity profile in the layer as 


u uty")... 9 U,V/ Vy (3.13) 


Substituting Eq. (3.13) into Eq. (3.12) and integrat 
ing it with respect to y, the deflection angle at the outer 
edge of the inner viscous layer 6, is 


0, —(ryye*/u,) (u,'/U,) (3.14) 


where prime denotes the differentiation with respect 
tox. Also, substituting Eq. (3.13) into Eq. (3.11) and 
integrating it with respect to y, the following equation 
holds: 


fPdy* — —(y,* pet) (dp, dx) i 


(u,/ Vy) [(7, Te) — 1] (3.15) 


By differentiating Eqs. (3.13) and (3.9) with respect 
to x at the outer edge of the layer, we have an equation 


(1/ py U,?) (dp, dx) — N(l i vrY.*) (u-?/U,*)0, ( U, 


ahs 
Fy (1 — a)" | [fd Y*/(f? — K)*/?] — [2y/(y — 1)] [mo/(1 + mo)] (1 — A) (U,2/ U0?) X 


>— K)'?] + My + 1)/(y — 1)) [mo/C + mo) }? + [C2 — 129/01 + mm)?! X 3.10 


K [(1 + mo)/mo] (Uo? U,*)A 


to describe the growth of y,*, 


(ue'/Ue) = (u,’/u,) + (fe’/fe) (3.16 


and an equation corresponding to Bernoulli’s equation 


(ue’ te) + (1/pyu,?) (dp./dx) = 
(Pe/ po) (Ue/tte)? (fe! /fe) (3.17 


where L’, is the value of U’ at Y* vo" 


Flow Conditions at the Edge 
Here a question arises 

Eq. (3.15). 

viscid layer. 


Consider the vorticity 2 in the outer in- 
The following relation* holds there: 


Q/p = function of ¥ alone 


Then, hypothesizing that the shear stress is proportional 
to the vorticity and the proportional constant is in- 
variable along a streamline, and assuming the constant 
shear stress in the layer of the undisturbed flow state 
we set up 


Te/ Pw = (Pe/ po) U,? (3.18 
Substituting Eqs. (3.14) and (3.18) into Eq. (3.15), we 
have 


Ve ¥,*) (ue U,) [(pe po) — (ue2/U,*)] (1 fe) 


Y,* (3.19 
where N (1/ Ye) [ fd Y* 
by virtue of y,* Y,* whereby u,/U, =u, U,. From Eqs. (3.14), (3.16), and (3.17), 
(1/pyU,2) (dp./dx) — (U,/tw¥.*) (ue?/U2)0. = [(pe po) 7/7 — (u.2/U.2)] (fe! /f,) (3.20 
Introducing A in Eq. (3.5) and A in Eqs. (3.10), Eqs. (3.19) and (3.20) are 
(1/2) (1 — A)7'49-Y) (dn K /dx) — [mo/(1 + mo)] (U2/Us2) (f2.W1 — A) (U,/ ae ¥-*) X 
N [1 — (K/f.2)]3/? (0,./A) =[(1 — A)78/7O-PI V1 — A] (U, oe ¥e*) {(1/f2) + [1/(y - 1)] X 
[mt (1 + mo)] (U2/U0)} [1 — (K/f,2)]¥2 (3.21 


where we use the approximate relation 


(pe/Po) — (u.?2/U.2) ~ [11 + mo) /mo| (Uo?/ U,2) [A/(1 — A) 


and 


(1/2) (AQ — A)~'9-YI (din K/dx) — [mo/( + m)] (UL 


From Eqs. (3.4) and (3.10), we have 


[m1o/(1 + mo)] (U.2/Uo?)t (3.22) 


Us?) (f2/W1 — A) (U,/vw¥.*) X 


[1 — (K/f.2)]8!2 (0,./A) = dInf,/dx (3.23) 


that is, the values of 7, in 
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g,/A (V Mo? — 1/2mo) — [(1 + mo) /2mo] (9, /Uo) (Up? /U,) (A — A428 -P40-DI Ryd in K/dx) (3.24) 
by using the approximate relation py ‘p, 1 + mo. 
Eliminating dx from Eqs. (3.21), (3.23), and (3.24), we have 
(1 — A) (d In f,2/d In K) D,/Ds 
where 
pD ‘1 + (1 — A)~° D Fst $1 + [1/(y — 1)] [mo/Q + mo) (U2 Ue) f 2 - 
, 9 l , ‘ a 7 5.25) 
(1/2) [V Mo? — 1/(1 + mo) (U.2/U,?) (1 — A) 1 — N)fA(f.’ K) 
D + [1 /(y — 1)] [nr /A + mo)] (U2 Uy) f2 + (1/2) (WM? — 1/0 + m)) (L Uy*) X 
(1 — A)’ -) NF2 (f.2 — K) 
ind F3 ((f.2 — K)*/2/V.*f.| Fs 
[he initial value of Y,* is decided by the condition as 
D; Oat K = 0. Next, introducing the nondimen- 
sional length val a 
g* = U.x/s (3.26) a 
: * aii 
we have the following equation to describe the pressure Po oat 
distribution at the wall: 
14} 
dx* l l fog 3 dD; | 
d i ? — f.- — ie ) 9 oO") 5 
Vin & VI-A Vi KM Ds (3.27) FAGE & SARGENT’S 
r 2 . Ee J 
where D; = 1+ N(1 — A)?7"F, 12) ne 
— THEORY 
Comparison With Experiments 
Compare the pressure distribution at the wall ob- a 
tained by the theory with the one measured by Gadd! 10 0.0 10 20 3.0, 40 
ry . U+ I . 
when JJ) ~ 3.0, the Reynolds Number depending on —— se” 
the boundary-layer thickness in the undisturbed flow Fic. 9. Pressure distribution at the wall when 1) = 1.472 and 
state do, Rs x 2.3 X 10 
R;, = Uobo/% = 7.5 X 104 
and a shock wave generated by a wedge in the free a 
boundary-layer assumptions. The same method is 


Assuming the re- 
120 
from Fig. 5 in the theory, a comparison between them 
Due to the flow separation the 


stream impinges on the flat plate. 
covery factor = 0.88, w = 0.76, and U,?/U,? = 


is shown in Fig. 8. 
pressure rise suddenly becomes milder after the vicinity 
of py» po ~ 2.0 in Gadd’s experimental data, as he has 
pointed out. 

Next, compare the pressure distribution at the wind- 
tunnel wall when the normal shock wave generated in 
the induced wind tunnel interacts with the boundary 
layer on the tunnel wall in the case of Jp 1.472, 
measured by Fage and Sargent.°® In the theory, we take 
R;~ 2.3 X 104, then lL? U2 ~ 480, by assuming that 
the stagnation states in this case are the ones of the 
standard atmosphere. Fig. 9 shows the comparison be- 
tween them where the origin is positioned at the point 
called ‘‘toe’”’ in their paper. The theory gives some 
larger values than the experiment with the progress of 
«*, perhaps due to the occurrence of flow separation in 
reality. However, they agree fairly well. 


CONCLUSIONS 
We describe the development of the laminar boundary 
layer upstream of the foot of shock wave by using the 
momentum and energy-integral equations on the usual 


applied to the flow downstream of the point of the 
boundary-layer separation, introducing a flow model 
as a jet-like flow. The results obtained in the theory 
agree fairly well with the experimental data as shown 
in Figs. 2 and 3. 

However, the method of approaching the subject 
in the case of turbulent flow along the similar lines as 

above theoretical on 
that is, the solution on the usual boundary 


mentioned has a contradiction 
the way 
layer assumptions throughout the layer has not always 
the same meanings in reality as shown in Fig. 7. 

In such circumstances, the turbulent boundary layer 
under consideration is divided into two regions, one the 
inner viscous layer and the other the outer inviscid 
layer. Then we determine the growth of the inner 
viscous layer, where the velocity profile is assumed 
f(y*), by the flow conditions at the inner edge 
The results obtained thus 


u/u, = 
of the outer inviscid layer. 
have fair agreement with those of the experiments as 


shown in Figs. 8 and 9. 
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Some Panel-Flutter Studies Using Piston 
Theory 


DAVID JOHN JOHNS* 
Sir IW. G. drmstrong Whitworth Aircraft, Limited 


SUMMARY 


The use of piston theory was recently advocated! for super 
sonic aeroelastic analyses, including the problem of panel flutter, 
ind this has stimulated the investigation reported here 

Linear piston theory is mainly considered, but some effects of 
introducing higher order terms are discussed 

Flutter of rectangular simply supported panels and of ellipti 
ily shaped clamped-edge panels is considered, and some justi 
fication is provided for the use of ‘‘static’’ aerodynamic forces 
ind the neglect of aerodynamic damping. Hence, it is concluded 
that Ackeret loading gives more exact results than piston theory 

Solution of the flutter equations is made by applying Galerkin’s 
method to a Ravyleigh-tvpe analysis using assumed modes of 


deformation 


SYMBOLS 


= speed of sound in undisturbed medium 
= modal coefficients defined in Eqs. (6) and (24 
= semiaxis of elliptic plate in chordwise direction (a 
= semiaxis of elliptic plate in spanwise direction (4 
= panel chord ( =2a for elliptic plate ) 
D = local flexural stiffness, Ef3/12(1 — v? 
E = Youngs modulus of elasticity of panel material 
= tension parameter, F/mc*w,* 
F = external force/unit width acting in mid-plane of 
panel (tensile force positive) in chordwise direc 


tion 

= structural damping coefficient 
Gn = matrix elements defined in Eqs. (15a) and (26a 
h = thickness of panel 
H,, H. = coefficients in assumed modes for elliptic plate 
k = reduced frequency, cw/2u 
k = stiffness parameter, cw /2u 
m = panel mass per unit surface area 
M = Mach Number 
N = coefficients in assumed modes for rectangular panels 
Pix = net perturbation pressure, positive downward 
Pniz = pressure coefficients associated with mode shape Z, 
p = pressure in undisturbed supersonic stream 
p = pressure on lower panel surface 

= dynamic pressure (1/2)pl? 
oo = coefficients defined after Eq. (28 

= time 
l = velocity of supersonic stream 
V2, V; = factors occurring after Eq. (32) 


= velocity of air normal to surface 
x, ¥, 2 = coordinates defined in Figs. 1 and 2 
Zi = vertical displacement of panel 
Zi = flutter mode shape 
Z = nth natural mode shape for rectangular panel vibrat- 


ing in vacuo 


The author would like to express his appreciation to the Direc- 
tors of Sir W. G. Armstrong Whitworth for permission to publish 
this paper 
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a = initial incidence of rigid panel 

8 = Vv i? l 

o, 6 = functions of /, k defined after Eq. (13 

im = panel-air mass ratio 

v = Poisson's ratio 

p = density of undisturbed supersonic stream 

a = density of pane! material 

w = frequency of oscillation 

w) = frequency associated with mode shape Z 

Q = frequency ratio squared (w;/w)*; in flutter calcula 
tions (w:/w)? (1 + ig 


= specific-heat ratio of air 
T = ratio of length of circumference of circle to its 


diameter 


(1) INTRODUCTION 


Feuer AND ZARTARIAN! have shown that simplifi- 
cations in the solutions of high-speed unsteady 
aeroelastic problems can be achieved by the use of 
“piston theory.’’ This term refers to any method for 
calculating the aerodynamic forces on a surface by ob 
taining a point-function relationship between the local 
pressure generated on the surface and the local normal 
component of fluid velocity at the surface, in the same 
way that these quantities are related at the face oi a 
piston moving in a one-dimensional channel. Essen- 
tially, the forces so calculated are ‘‘static forces.” 

In reference 1 the use of piston theory has been 
extensively explored and its limitations defined, and it 
is with this background of information that the present 
investigation was started. 

The panel-flutter problem has received considerable 
treatment in recent years, and it is interesting to note 
that the majority of the investigations have considered 
unsteady flow aerodynamic theory rather than so-called 
‘static’ forces. (For an extensive list of references see 
Hedgepeth’s recent paper.”) 

Typical of such studies is that by Nelson and Cun- 
ningham* in which they obtained the pressure on the 
upper surface of the panel by the theory for linearized 
unsteady supersonic flow and introduced perturbation 
pressures on the lower surface from acoustical theory. 
Their results showed good correlation with experiment 
for the low Mach Number region (J < 2.0). 

Hedgepeth,” however, has assumed “‘static’’ air forces 
and neglected aerodynamic damping effects, justifying 
this on the grounds of experience. He considered both 
aerodynamic strip theory and aerodynamic surface 
theory and obtained excellent agreement between the 
two analyses. For strip theory, the lateral loading is 
given by the simple Ackeret value, to which linear piston 
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Fic. 1. Coordinate system of flexible two-dimensional panel. 


theory reduces, for large .\/ and neglecting aerodynamic 
damping. 

In this present paper, therefore, we again consider 
the problem of the flat unbuckled panel and, using 
linear piston theory, derive the flutter boundaries for 
rectangular simply supported panels and clamped-edge 
elliptic panels, for Mach Numbers greater than 7/2. 
In the former analysis, the justification for neglecting 
aerodynamic damping forces is obtained, and, as would 
be expected, the results then compare with Hedgepeth.’ 


(2) STATEMENT OF THE PROBLEM 


The configurations to be considered here are shown 
in Fig. 1 for the rectangular simply supported panel, 
and in Fig. 2 for the elliptic clamped-edge panel. A 
flat plate of uniform thickness / is mounted in a rigid 
wall with air flowing over the top at a supersonic Mach 
Number 1/7. The plate is subjected to midplane stresses 
and is assumed to be undergoing simple harmonic 
motion of frequency w. 

In the following analysis, small-deflection thin- 
plate theory is assumed. The differential equation of 
motion for any such plate in general may then be writ- 


ten as 


DV zi + m(0?s/Ot?) — F(0?s/0x?) + 
(bp. —fP) + Pen =0 (I) 


where the vertical displacement Z,, », may be expressed 
asZ,.e™. 

The boundary conditions applicable will, of course, 
differ for the two cases to be considered. In this anal- 
ysis the constant pressure term (Pp. — p;) of Eq. (1) 
is considered to be zero. This term, if included in the 
analysis, would only add a particular solution to the 
results, and hence its exclusion in no way affects the 
generality of the results for the oscillating plate. 

Using piston theory,' the exact expression for the 
instantaneous pressure p on the face of the piston is 


P/ha = {1+ [(y — 1)/2] (w/a) }?”” 2) 
which may be approximated by the linear relation 
Pp — pa = paw (3) 


or the second-order expression 











SCIENCES~NOVEMBER, 1958 
= in subs 
Zz, 4y ind app 
yatiol 
nd the 
minant 
Fic. 2. Coordinate system of flexible elliptic panel The1 
oe i 29 (en /7 | 
p Pp pa (w/a) + where 
[(y + 1)/4] (w/a.)*} (4 
The velocity w at any element of a panel whose modal 
shape is Z,,) is given as 
wW U(0z/0x) + (O2/d0t) + Ua 5 
where the term Ua is a measure of the initial incidence 
of the rigid panel. It can be shown that this term may 
be neglected, if linear piston theory is considered, as it 
introduces another constant-pressure term only. Its | and wh 
inclusion in the second-order expression of Eq. (4) intro- 
duces further time-dependent pressure terms, the Z 
effects of which will be discussed later. st 12 
With hi 


Actually piston theory is really valid only for Mach | 
Numbers greater than 1/ 2.5 (see reference 1). Its | 
use will be considered in this note for .1/ / 2, 2, and 
6 because of the simplicity with which aerodynamic 
damping forces can be included, and hence for all these | where 
Mach Numbers it will be shown that such forces may 
be neglected. 
The 
(3) ANALYSIS 

(3.1) Rectangular Simply Supported Panels 

A Rayleigh-type analysis—similar to that of refer- 
ence 3—1is made, the same notation being adopted for 
convenience, with the coordinate « nondimensionalized 


ve — , (1.1) 
rhe flutter mode shape Z,,) is approximated by a 


linear combination of the form It 1s 
. compl 
ea 2 eS ae (@) where 
n=1 efficies 
where the coefficients a, may represent complex ampli- | ‘ 2” 
tudes and the functions Z,:,) are the mode shapes for | Stabili 
the panel vibrating in a vacuum with no axial force F | /: k,a 
acting. The shapes Z, satisfy the pinned edge bound- only a 
ary conditions and the differential equation soluti 
- quadr 
~4\ 7 2 y — 7 ° 
(D/c\Z, ° — wx*mZ, = 0 (7) soluti 
For a two-dimensional panel, Eq. (1) may be rewritten 
. (3.1.2) 
as 
” Fro 
(D/c)Zq* — mw’Zy) — far 
= 1S ¢ 


(F/e*)Za)"' + po, ne” =0 (8 





PANEL-FLUTTER STUDIES USING PISTON THEORY 6S 
in substituting Eq. (6) into Eq. (8), utilizing Eq. (7), TABLE | 
7? nd apply ing the Galerkin process,’ we obtain .V linear Mode Shape Factors, Eigenvalues, Frequency Ratios, and 
; a Stuctural Integrals for First Four Normal Modes 
equations 
GuGe . . Gin in \°| Motes Ni Ke wali Aum B | D 
Gay Go» . . Go ade = 0 1 | - l 0.5 0.5 0 1/3 0 8/1 
(9) 2 ] 2 j 0.5 2.0 1/3 0 12/5 0 
. | 1 Q 0.5 1.5 0 12/'5 0 24/7 
Gy; Gyo ; ; Gwyx ay QO 4 ] 4 lf 0.5 8.0 8/1 0 24/7 0 
“4: : Nort 1 B 0 ~ 
nd the flutter condition may be expressed in the deter- 
minantal form 
Ginn| = 0 (10) 
| 
4 Y y I ] 
The matrix elements are given by ? L/4.Mk*, D, i, Linley AN 
HtAmn — Q[(wn?/or?)Amn + fBmnly — Cnn (11) § = 1/2Mk, k cw/2U 
a where The solution of the flutter determinant Eq. (10) can 
{ t a . . - 9 
' ae now follow established lines (e.g., see reference 3). As 
modal sii oe saanate in reference 3, no more than four modes will be con 
, sidered, the mode shapes Z,, and their associated natural 
: Q = (w/w)? B [ Z,17 Jdx frequencies w, being obtained from Eq. (7). 
3) mr m ? ° 
—— For the pinned edge plate 
dence >] 
1 may f = F/mctw,*, CL. = | ZmPr(x)dx Z; Nisin K,x | 4 
ms J0 P f (14) 
, AS Wn K,.2 V D/mc*\ 
Its | and where 
ntro- The values of the various parameters and the mode 
r . a tut i —_ tw 2. 3 E ‘i . - ° 
the Zi, Zine”, Pn(X) = Pe, of on shape integrals for the first four normal modes are given 
rea ae . in Table 1. 
With linear piston theory ” I ‘ ale : ; : 
Mach Substitution of these values into Eq. (10) leads to 
Its Pu. t) = (pu, M) [U(0z/0x) + (02/04) ] (12) the flutter determinant shown below: 
_ and 
amic Ca , = oD, oT 6iA,, n (13) Gy Gy Gi Gi 
hese where Gai C22 G23 C74 0) (15) 
mav G31 G32 Ge3 G4 
Guy Gas G4: Gs 
The terms in the flutter determinant can be defined as follows: 
Gy = 0.5[{1 — QC + rf) — (26/p)], Gy = —Gay = (4/3) (o/pn) 
Gee 0.5[1 — Q(116 + 42?f) — (16/p)], Gog = —Gg = (12/5) (b/p) 
hae. G33 = 0.5[1 — Q(81 + 92?f) — (26/pn)], Gu —G; (S/15) (@ pw) (15a) 
7 Gu = 0.5[{1 — Q(256 + 16r°f) — (76/u)], Gay = —Gig = (24/7) (6/u) 
zed Gis G31 = Go = Ge 0 
(3.1.1) Solution of the Flutter Determinant 
ya 
It is convenient to interpret the 2 of Eq. (15a) as the considered : 
complex quantity (w;/w)*? (1 + 7g) rather than (#)/w)* (1 + wf) (16 + 4r2f)Q2 — (17 + Sef) [1 — (6/u) 2+ 
(6 where g may be regarded as a structural damping co- {11 — (i6/u)|2 + (64/9) (¢?/n?)} 0 (16) 
efficient. To solve the determinant as shown, recourse 
. ‘ : ‘ P we ace , ‘ »O = », /@ 2 +. 4 ; i 
ji- | to an electronic computer is desirable, and hence the If we assume that one 2 (w/@)* (1 8) solution 
or | Stability boundaries can be obtained for ranges of .1/, of Eq. (16), for any value of structural damping coeffi 
F | fk, and g (cf. reference 3) If. however. we consider cient g, may be made. The following stability boundary 
d- | only a two-mode solution (in reference 3 only two-mode is obtained for g = 0: 
solutions have been made for pinned edge panels), the 7 5+ 4 TP 36M! ™ 
quadratic equati ained is amenable analytic: =i. er (1¢) 
7) ratic equation obtained is amenable to analytical ‘ 17 + 5rf (16/9) — R? 
solution. 
“1 and the frequency ratio 
(3.1.2) Two-Mode Analysis P 
®1/0 = vz (17 + 5r°f) (1S) 


From Eq. (15) the following quadratic equation in 


2 is obtained if the first two normal modes only are Using Eqs. (17) and (18), the flutter boundary de 
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STIFFNESS PARAMETER, 2k, 
Fic. 3. Stability boundaries for simply supported, two-dimen- 
sional panels (f = g = 0, 2 mode, ---- - 4 mode). 


fined by the variation of | /u with the stiffness parameter 
2k; can be obtained. Such boundaries are shown in 
Fig. 3 for f = 0. It is worth noting that it is necessary 
to solve Eq. (17) only once for any particular value of 
k, the corresponding values of 1/4 varying linearly with 
M. 

Also shown in Fig. 3 are the hyperbolas (1/u) X 
2k,, for aluminum panels at sea level and at 50,000 ft. 
for the Mach Numbers y/2, 2, and 6. The product 


(1/u)-2k: = [K2/V24(1 — »2)] x | 
(p/a) V pE/ ag (19) 
= (0.01385 for aluminum panels at 

M = V2 sea level 


The intersection of the hyperbolas with the corre- 
sponding stability boundaries in Fig. 3 gives the critical 
values of 1/u = (p/oc)-(c/h) and, hence, of h/c, for the 
different flight conditions. 

Inspection of Fig. 3 shows that even for aluminum 
panels at MJ = +/2, sea level, the maximum value of 
2k; is 0.16 which corresponds to a maximum value for 
k of 0.234. It is seen, therefore, that exclusion of the 
k* term from Eq. (17) would introduce a maximum 
error of less than 2 per cent in the value for 1/u. Ata 
Mach Number of 6 where piston theory is more valid, 
the error is considerably smaller. 

For other materials the corresponding errors will 
depend on the value of the parameter VE/o8 [see Eq. 
(19) ]. 

It should now be noted that exclusion of the k® term 
from Eq. (17) is equivalent to neglecting aerodynamic 


SPACE SCIENCES 


NOVEMBER, 1958 


damping. Hence the smallness of the error introduce, 
justifies the assumption of Hedgepeth* that dampin 


forces be ignored. Eq. (17) then reduces to 


1/uk? = [((5 + wf)/(17 + 5xf)] (9102) (99 
from Eq. (18) 

ki/k =V2/(17 + 5r*f) 21 
therefore, 

1/uky? = (94/4) (5 + wf) (22 


and this equation can be rearranged to give 

2qc*/MD = (974/16) (5 + xf) (23 
Eq. (23) agrees with the result obtained by Hedgepeth 
with I substituted for 8 = VM? — 1. 

It should be noted that the above results apply only 
to the buckling of the plate. Neglecting the effects oj 
aerodynamic loading, this corresponds to a value oj 
m*f = —1 for a two-dimensional simply supported plate 

To summarize the results of this analysis, Fig. 4 gives 
the variation of critical thickness ratio hc with Mach 
Number for unstressed aluminum panels using Eq. (23 
For comparison, the results of references 2, 3, 4, and 5 
are also plotted. The correlation between the various 


results is excellent. However, since aerodynamic 
damping can be neglected, Ackeret loading is more 


exact than piston theory. 


(3.1.3) Four-Mode Analysis 


To solve the full quartic equation obtained from 
the flutter determinant Eq. (15), recourse must be made 
to an electronic computer, in which case, the deter- 
minant may be solved in the way described in reference 
3. 

A thorough analysis of the quartic has not been 
deemed necessary, but enough solutions have been 
made to confirm the results obtained from the two- 
mode analysis and to show the effects of increased 
number of modes. 

The quartic was solved for one Mach Number only, 
MT = 2, and the stability boundary obtained is shown 


dotted in Fig. 3. It is seen that increasing the number 
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PANEL-FLUTTER STUDIES 


of modes from two to four has moved the stability 
boundary to the left-—1.e., the two-mode solution gives 
a conservative estimate for the critical panel thickness. 
The effect of this conservatism on the value of critical 
panel thickness is shown in Fig. 4 by the results of this 
analysis and those of reference 2. Similar effects were 


also obtained from the analysis of reference 3. 
3.2) Elliptic Clamped-Edge Panels 


Up to the present time, most panel flutter analyses 
have dealt with rectangular panels with various edge 
restraints. For design purposes, any particular panel 
of arbitrary shape has had to be approximated to an 
equivalent rectangular shape, with the appropriate 
choice of boundary restraint. 

Because panels are often circular or at least elliptic, 
and to help in the design of panels of arbitrary shape, 
the investigation reported below has been made. 

A Rayleigh-type analysis similar to that in Section 
3.1) has been made. 

The flutter mode shape Z,,) is approximated by 


Zic. t) = [a(Zi1/M) + a2(Z2/ He) je (24) 
where 


Z, = [(x?/ae) + (y?/b?) — 1}? 
Z2 = x[(x?/a?) + (y?/b?) — 1]? 
H, = [(3/a*) + (3/64) + (2/a7b?)] 
H, = {(5/a*) + (1/b*) + (2/a7b?) | 


and where the coefficients a), @» may represent complex 
amplitudes. 

Actually it is more usual to assume, as modes of dis- 
placement, the mode shapes for the unstressed panel 
vibrating in a vacuum, but these are unknown to the 
author. The modes assumed have been obtained by 
Bryan’ and are respectively due to pressure loadings of 
the form p (const.) and px. 

Therefore, assuming that mid-plane stresses are zero, 
the following equation of motion for the elliptic panel is 
obtained. 

DV4z + m(0°2/O0t?) + pa,» = 0 (25) 
and, hence, substitution of Eqs. (24) and (12) into 
Eq. (25) yields, after applying Galerkin’s method, the 
resultant flutter determinant. 

Gy Gis 


Go, Go 


=e §) (26) 


The limits of integration for the various structural 
integrals are 
~b V1 — (x?/a2) < y<t+b V1 — (x? a’) 
—a<x<+a 
and where 
Gy = (8D/3) + [(pUie — Mmw*)/5MHA;]| 
Gy» = a2D + a*[(pUiw — Mmw?)/60M 2) 


Gy» = pU?/10MH 


(26a) 
Gx = —(pU?/10MAj) | 


Solution of Eq. (26) for the case of zero structural 
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damping yields the instability boundary 
1/u? = [1 — (r/75S?)] [4.17°R* (3 — R*)] 27) 
and the flutter frequency 
w = VI5SD pap 28) 


where ¢ Si\f12./3 and S (2/7, + 9s) 45 are 
geometrical functions of the plate. 

If aerodynamic damping is neglected, and this may 
be justified as in Section (3.1.2), Eqs. 27 and 28 reduce 
finally to 
(SO +1/3/9) X 


[39 + 14(a? b*) + 3(a* b*)] (29 


296 , MD 


As pointed out for the rectangular plate, it is probably 
more correct to replace / by V i? — 1, thus using 
Ackeret loading. 
From Eq. (29), the critical thickness ratio for panels 
of aspect ratio | and infinity are, respectively, 
a/jb=i1: h/e 0.295 ¥ qd MEI 


a/b =0: h/« 


(30) 
0.331 Vq/ME\ 


The ratio of these results, 0.351 0.295 = 1.128, agrees 
well with that obtained by Hedgepeth? for rectangular 
simply supported panels. This shows that the critical 
panel thickness ratio for unstressed panels does not vary 
much as the panel aspect ratio increases above unity 
for Mach Numbers greater than about v 2. 

The result for the two-dimensional clamped edge 
panel—.e., 2gc* JD = 600—compares with the pre 
vious result obtained in Eq. (23) for the simply sup 
ported panel—-i.e., 2gc* WD 274. Therefore, the 
ratio of critical thickness ratios for simply supported 
and clamped edge two-dimensional panels is W600 274 
= 1.3. This ratio agrees approximately with that ob 
tained in reference 3. 

From Eq. (30), the critical thickness ratios for two 
dimensional aluminum panels for 1J = +72, at sea 
level and 25,000 ft., are, respectively, 0.00375 and 
0.00269. These results compare with values ot 0.0045 
and 0.0030 obtained by Nelson and Cunningham? using 
unsteady flow aerodynamics. If, however, Eq. (30) is 
modified to include the term V2? — 1 rather than V/, 
the values obtained for critical / c are 0.00421 and 
0.0030. Hence, it is concluded that Ackeret loading is 
more exact than piston theory, especially at low Mach 
Numbers for panel flutter analyses where aerodynamic 
damping may be neglected. 

Therefore, using Eq. (29), Fig. 5 presents the vari- 
ation of the parameter 2gc* 17D with panel eccentricity 
a/b. It can obviously be concluded that as a/b > « 
(i.e., a > 0) then h/c — 0, and the most critical value 
of h/c is for the two-dimensional panel. This is also 
shown by Fig. 6 which gives the variation of the critical 
panel thickness ratio with Mach Number for the two- 
dimensional (a/b = 0) and circular (a 6 = 1) clamped 
edged panels. Fig. 6 has been derived from Eq. (30) 


with / replaced by VM? — 1. 
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Fic. 5. 


(3.3) Effects Due to Initial Incidence of Panel 


The second-order expression for the pressure at a 
point derived from piston theory is given in Eq. (4) 


bp — pe = pa.”\(w/a.) + [(y + 1)/4] (w/a.)?} 
where w is given by Eq. (5) 


w = U(dz/dx) + (02/dt) + Ua 


VI1Z., 


Because the U, term is not time-dependent, its inclu- 
sion with linear piston theory only adds a particular 
solution to that of the differential equation of motion 
for the oscillating plate and may therefore be neglected. 

If, however, second-order piston theory is considered, 
the U’, term introduces extra factors into the flutter 
equation. Substitution of Eq. (5) into Eq. (4) pro- 
duces the following expression for the pressure in the 
flutter equation: 


[Pur, ole = (0U/M) [U(Oz/0x) + (02/0t)] X 
) 


11 + [Maly + 1)/2]} (31) 

















i.€., [Pre, » Je = [Pr, V2 (32) 
sr T 7 , 
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where [/:,, » |:,2 are pressures due to linear and second 
order piston theory and 


Vs = ;1 a [Maly + 1) 2k 


Introduction of the factor V into the flutter equations 
produces the result that the new stability boundary 


[1/ulo = (1/V2) [1 / uh (33 


where [l/u]:,2 are the critical panel-air-mass ratios 
using linear and second-order piston theory. 
Similarly, if the third-order expansion of Eq. (2) js 


obtained—viz., 


p— pa = pa."\(w/az) + [(y + 1) 4] X 
2+ [(y + 1)/12] (wa.)*} (34 
the factor 


V3 = 1 + [Waly + 1) yA + 
[. 


WPat(y + 1)/4]} 39 


Inspection of Eq. (33) shows that the factor V2 can 
be stabilizing or destabilizing, depending on the sign 
of a. It should be remembered, however, that there 
is a physical limitation in the use of piston theory in 
that the surface involved should not be inclined too 
sharply to the direction of the free stream. This limi- 
tation implies that .\Ja < 1 (see reference 1), and, hence 
it would seem that l’; = | for all practical purposes 
The conclusion can therefore be drawn that initial 
incidence of the panel has no first-order effect on the 


stability of the panel. 


(4) CONCLUSIONS 


The use of piston theory in panel flutter analyses 
has been considered, and the flutter boundaries for 
flat, rectangular, simply supported panels and clamped- 
edge elliptic panels have been obtained. Justification 
has been made for the use of ‘“‘static’’ air forces only, 
neglecting aerodynamic damping, and it has been con- 
cluded that Ackeret theory is then more accurate than 
piston theory, especially at low Mach Numbers. 

The results for the rectangular panels agree well with 
those of other studies using more complex unsteady 
aerodynamic forces. 

The results for the elliptical plate are the first known 
to the author for other than rectangular panels. Al- 
though a two-mode analysis only was made, the accu- 
racy of the results would appear reasonable if the limit- 
ing case of the two-dimensional panel is compared with 
more detailed studies. 

In the above investigations linear piston theory was 


assumed. However, an analysis of the effects of initial 
incidence of the panel using higher-order piston theory 
terms would appear to indicate that the flutter bound- 


ary is not a function of the incidence of the panel. 
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Magnetohydrodynamic-Hypersonic Flow Past 


a Blunt Body 


WILLIAM B. BUSH* 
The Ramo-Wooldridge Corporation 


SUMMARY 


The flow field near the stagnation point of an axisymmetric 


blunt body in hypersonic flow is analyzed for a spherical detached 
shock 


which is assumed to be incompressible, inviscid, and of constant 


\ magnetic field is applied to the fluid in the shock layer 


electrical conductivity. 

\ family of solutions is found for which the body shape is a 
sphere concentric with the detached shock. For this family of 
solutions, the standoff distance and pressure relief increase with 
increasing magnetic field strength. Moreover, with increasing 
magnetic field strength, the gradient of the tangential velocity 


it the body decreases. 


SYMBOLS 


B = magnetic field 

B,, Bo = magnetic field components in 7, @ directions 

( = radius of spherical shock 

E = electrical field 

e,,@9,@4 = unit directional vectors in r-, 6-, ¢-directions, re- 
spectively 

f(r = function of radial distance introduced by the ve- 
locity vector 

F(R) = f(r)/Uc? 

J = electric current density 

m(? = function of radial distance introduced by the mag- 
netic field 

W(R) = m(r)/ac* 

p = pressure 

r = nondimensional parameter associated with the 
pressure on the body 

Vv = ga’c/p.U 

q = velocity vector 

u, 1 = velocity components in 7, 6 directions 

l = velocity in the free stream, parallel] to axis of sym- 
metry 

V, = nondimensional velocity gradient at the body 

r, 0, @ = coordinate axes 

R = 7/C 

Rey = magnetic Reynolds Number, oucl 

a = constant associated with the magnitude of the 
magnetic field 

A = nondimensional standoff distance 

€ = Po./Pp 

¢ = component of vorticity in ¢-direction 

m = magnetic permeability 

p = density 

a = electrical conductivity 

® = meridian radius 

Q = vorticity 


= quantity in the free stream 
= quantity right behind the shock 
= quantity normal to the shock 

= quantity tangential to the shock 


= quantity evaluated at the body 
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) = stagnation quantity in free stream 
= jump in quantity through shock wave, 


INTRODUCTION 


HIS PAPER? considers the steady-state problem of 

hypersonic flow over an axisymmetric blunt body 
(Fig. 1). Due to the bluntness of the body, there is a 
detached shock, and attention is focused on the gas 
that lies in the shock layer between the detached shock 
and the nose of the blunt body. 

In the free stream the fluid is considered to be in- 
viscid and to have a velocity, l’, uniform and parallel 
to the axis of symmetry. Due to the hypersonic 
speed of the flow in the free stream, the Mach Number 
in the shock layer near the nose of the blunt body is 
quite low, and the approximation is made that in this 
region the fluid behaves like an incompressible fluid. 

Although the electrical conductivity of the fluid in 
the free stream may be considered to be negligible, the 
temperature rise across the shock wave is large enough 
to produce a degree of ionization in the shock layer, 
which means that the gas here should be considered an 
electrically conducting fluid and capable of being in- 
fluenced by a magnetic field. Since this is an intro- 
ductory paper, the problem is simplified by making 
the assumption that the electrical conductivity of the 


+ For an analysis similar to the one carried out in this paper 
and a discussion of earlier work done on the stagnation point 
problem with a magnetic field applied, the reader is referred to 
the note of Kemp,! which has appeared while the present paper 


was in the process of review and revision. 





uv, 





Fic. 1. Diagram of blunt body and detached shock. 
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gas in the shock layer is constant. 


sumptions that the heat conduction and viscosity effects 


may be neglected in the shock layer are also made. 


In this analysis, the shock front is spherical, and a 
body shape is found which is consistent with such a 
shock shape and with the flow conditions in the shock 


layer. 


specified in the free stream and the magnetic field found 


at the body. 
sidered is the one due to a dipole at the origin. 
means that the electrical coils in the body are arranged 
in such a manner that the applied field, modified by the 
flow in the shock layer, appears to the free stream as 
a field due to a dipole at the origin. 

Furthermore, it should also be pointed out that the 
resultant field in free stream is, in general, not of the 
same shape or magnitude as the field produced by the 
coils in the absence of the plasma currents. This paper 
deals primarily with the effect of the hydromagnetic 
interaction on the flow field and, thus, the means of 
producing such a magnetic field are not considered. 
Such a field satisfies the conditions that its divergence 
and curl be zero and that its magnitude approach zero 
as the radial distance approaches infinity. 

The electrodynamical and fluid-dynamical jump con- 
ditions across the known shock shape, with known up- 
stream conditions, make possible the explicit expres- 
sion of the boundary conditions just behind the shock. 


STEADY-STATE MAGNETOHYDRODYNAMIC 
EQUATIONS 


(1) THE 


The steady state fluid mechanical and electrody- 
namical relations for the flow field are expressed by the 
following equations:?~‘ 

V-(pq) = 0 (1.1) 


V(9g?/2) —qx(V xq) = 
— (1/p)V£ + (1/p) (jf x B) (1.2) 


VxE=0 (1.3) 
V-B= 0 (1.4) 
J=o0(E+qxB) (1.5) 
VxB=yj (1.6) 


These equations shall be considered for a spherical 
coordinate system where r and @ are the coordinates in 
the meridian plane and ¢ is the coordinate perpendicular 
to the meridian plane. 

If the electric current density has a component only 
in the ¢-direction, then Eqs. (1.2) and (1.6) show that 
the quantities B and g must lie in the meridian plane. 
This is consistent since the flow behind a spherical shock 
would be expected to have components only in the r- 
and 6-directions. 

Since B and g have components only in the meridian 
plane, their vector cross-product, g x B, has only a com- 
ponent perpendicular to the plane. Eq. (1.5) shows 
that E at the most, can have only a component in the 
¢-direction. Moreover, the curl of E must be zero. 


SPAC 


The further as- 


In the same manner, the magnetic field is 


The free-stream magnetic field con- 
This 
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VxE = (1 rsin@) [0(/, sin 6) /00Je, — 
(1/r) [O(Egr)/or]eg = 0 (17 


The two possible solutions are 


Fy, = (constant) (1,7 sin 6) and Es. = (12 


The second solution is the only satisfactory one, becaug 
of the singular behavior of the first solution wher 
6 is zero, leading to the conclusion that E is zero 


(2) THe BOUNDARY CONDITIONS 


Having expressed the equations that govern mag. 
netohydrodynamics, the next step is to specify th 
boundary conditions for the particular problem that js 
to be solved. 

A study of Maxwell’s equations shows that the mag 
netic field, normal to the shock at the shock front, is 


continuous. Thus, 


fi = 6 (2) 


If the current density parallel to the shock front does 
not become infinitely large, and if the magnetic perme. 
ability in the free stream is that of the fluid in the shock 
layer, then 

[B.] = 0 


(>> 


Since the conductivity of the fluid in the free stream 
is assumed to be zero, the curl of the magnetic field in 


the free stream must be zero. Thus, 
V-B=0 (2.3 
and Vz = 0 (2.4 


in the free stream. The other condition is that the 
field, B, vanishes for all @ as r approaches infinity 
Inspection of these equations shows that one possible 


solution is the dipole solution, 


(2ac* cos 6/r*) e, + (ac*® sin 6/r?) e, (2.5 


B= 


where a is a constant, having the dimensions of B, and ¢ 
is the radius of the shock. This is the behavior of B 
which is postulated in the paper. 

Making use of Eqs. (2.1), (2.2), and (2.5), 
(2.6 


B, = 2a cos Oe, + asin 6e, 


The fluid mechanical jump conditions across the 
shock are 


[pgn] = O (2.7a 
[a:] = O (2.7b 
[pgn? + p] = O (2.7¢ 


since there are no magnetic terms in these expressions 


because of the continuity of B across the shock. Across 
a spherical shock, Eqs. (2.7a) and (2.7b) yield 

q: = —«€U cos be, + U sin Oe, (2.5) 
where « = (p,/p). Eq. (2.7c) becomes 


ps = pj + p,U? cos? 6 (1 — e) (2.9) 
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Although it plays no part in the rest of the analysis, 
it is of interest to note that just behind the shock, 7, = 
—gal’ sin 6 cos 6 (1 + e) ¥ O, while in the free stream 
the current density is zero. 

Thus, the boundary conditions show what 6 depend- 
ence is required for B and g in the shock layer. 


(3) ANALYSIS 


The momentum equation for the nonviscous, elec- 
trically conducting fluid in the shock layer is 


V(q?/2 —qx(Vxqgq) = 
— (1/p)Vp + (1/p) (fx B) (3.1) 


Taking the curl of the momentum equation, consider- 
ing that the fluid is incompressible, one gets 


—V x (qx (Vxq)] = (1/p)Vx(jxB) (3.2) 
Introducing Q, the vorticity, defined by 
Q=Vxq = Oe, + 0e, + fe, (3.3) 
Eq. (3.2) becomes 
—V x(¢gxQ) = (1/p) Vx (jx B) (3.4) 


For ail axisymmetric, incompressible fluid, the left-hand 
side of the above equation may be re-expressed as 


[-—V x (gx Q)]-e, = @:Ve — (f/a)q:Vo (3.5) 


where @ is the meridian radius, the distance normal to 
the axis of symmetry in the meridian plane. In spheri- 
cal coordinates, @ = r sin 6. Therefore, the final ex- 


pression for the curl of the momentum is 


q:V¢ — (¢/a)q:Vo = (1/p)[V x (J x B)]-e, (3.6a) 


This is the first equation that is studied; the second 
equation comes from Maxwell’s equations and Ohm’s 
Law. Itis 


J = o(fE + qx B) = o(qxB) = (1/w)VxB (3.6b) 


Considering first the left-hand side, or fluid mechani- 
cal side of Eq. (3.6a), an expression for g must be ob- 
tained. Since the fluid is incompressible, the divergence 


of g in spherical coordinates is 
V-g = 0 = (1/r*) [0(r?u)/or] + 
(1/r sin 6) [O(v sin 6)/00}] (3.7) 
Rearranging, 
[O(r?u sin @)/Or|] + [O(r v sin 6)/06] = O (3.8) 
A solution that satisfies this equation and at the same 


time also satisfies the 6-dependence of the boundary 
conditions at the shock is 


q = [2f(r) cos 6/r*Je, — [f’(r) sin @/rje, (3.9) 


where f(r) is the function of r to be determined. 

Here it should be stated that, throughout the rest of 
the analysis, the terms of O(@)* shall be neglected in the 
expressions for sin 6 and cos @, both in the equations 
and in the boundary conditions—i.e., 


cos 6 = 1 — 6°/2 (3.10a) 
sin 6 ~ 6 (3.10b) 


It will be shown that with this approximation the 
boundary conditions are still satisfied, and, that, at 
the same time, the partial differential equations in 7 and 
6 can be reduced to ordinary differential equations in r. 

Retention of terms of higher order in 6 for sin 6 and 
cos 6 would, of course, satisfy the boundary conditions. 
However, because of the complexity of the equations 
for even this simplified problem, such retention would 
prevent the reduction to ordinary differential equa- 
tions. It should be noted at this point that a non- 
spherical shock front may coincide with the given one 
to the approximation cited. 

With this new limit on the problem, 


q = [2f(r) (1 — 6°/2)/r*] e, — [f'(r)-0/rjeg (3.11) 


Next, evaluating the vorticity expression, Q = V xq 
by introducing the expression for g in Eq. (3.11), it 
can be seen that 

¢ = —(6/r) (f’’ — 2f/r?) 3.12) 
where the primes denote derivatives with respect to r. 
The substantial derivative of the vorticity is 


q:Vé = u(O¢/Oor) + (z/r) (OF/08) (3.13a) 
= (6/r*) (—2f f'" + Off /r + ff" + 
2f f’/r? — 12f?/r®) (3.13b) 


Using Eqs. (3.11) and (3.12) and the definition of @ in 
spherical coordinates, the second term on the fluid 
dynamical side of Eq. (3.6a) is 


(¢(/@)q:Voa = —(@/r’) X 
(Qf f''/r —f'f"! — Af?/r® + 2f f'’/r*) (3.14) 


This means that the left-hand side of the equation 
for the curl of the momentum 1s 


[—V x (¢x Q)]-e, = —(2f/r*) X 


(f7? — 2f''/r — 2f'/r? + 8f/r*) 0 


vee 3.15) 

In order to consider the electrodynamical side of the 
Eq. (3.6a), there must be an expression for B. The 
field must satisfy the condition that its divergence be 
zero. In an analysis entirely similar to that for the 
divergence of the velocity, B can be expressed, as 


follows: 


B 


3.16a) 


[2m(r) cos 6/r?)e, — [m’(r) sin 0/r] e, 


v 


[2m(r) (1 — 6?/2)/r?] e, — 


[m’'(r)-6/r] eg (3.16b) 


where m/(r) is a function of 7 that is to be determined, 
but where the 6-dependence satisfies the boundary 
conditions at the shock. 

Since E is zero, Ohm’s Law becomes 


Jj = o(qx B) 3.17) 
and the curl of the Lorentz force is 


Vx(jxB) = oVx((gxB)xB 3.18) 
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Making use of Eqs. (3.11) and (3.16b), Eq. (8.18) be- 


comes 


{v x [qx B) x B}}-e, = 
(46m/r°) (—fm"’ + f’’m + 4fm’/r — 4f’m/r) (3.19) 


Therefore, Eq. (3.6a) becomes 


Sf" — 2f"/r — Af! /r? + 8f/r*) = (—2Zom/pr*) X 
(—fm'’ + f'’m + 4fm’/r — 4f’m/r) (3.20) 


Introducing the following nondimensional variables, 
R = r/c 
F(R) = f(r)/Ue? 


(3.21a) 
(3.21b) 
M(R) = m(r)/ac? (3.21c) 
where a was defined in the section on the boundary 


conditions, Eq. (3.20) becomes 
F(F’"’ — 2F"/R — 2F'/R? + 8F/R?) = 
— (2/R?) (ca®c/p,.U) (p/p) M (—FM”’ + 
F’'M + 4FM’'/R — 4F’M/R) (3.22) 
Defining the dimensionless parameters, 
Q = (ca’c/p_U) (3.23a) 
and = psp (3.23b) 
The curl of the momentum equation is 
F(F’" — 2F"/R — 2F’/R? + 8F/R*) = 
— (2/R?)eMQ(— FM" + FU'M + 
1FM’/R — 4F’M/R) (38.24) 
Turning to Eq. (3.6b), we see that Ohm’s Law may 
be expressed as 
a(qx B) = (1/n) (Vx B) (3.25) 
Substitution of Eqs. (3.11) and (3.16b) and the proper 
derivatives into Eq. (3.25) yields 
—(m"’ — 2m/r?) = (2ou/r?) (—fm’ + f’m) (3.26) 
By introducing the dimensionless quantities of Eq. 
(3.21), it can be seen that 


—(M" — 2M/R?) = (2/R?)Rey X 
(—FM’ + F’M) (3.27) 
where Rey = (oucU) 


Thus, there are two ordinary differential equations, 
Eqs. (3.24) and (3.27), for the two unknowns, F(R) and 
M(R). 

Now the boundary conditions must be re-expressed 
in terms of the new variables. First, considering the 
boundary conditions on the field, 


— @°/2)e, + a-be, = 
2aM(1) (1 — 6°/2)e, — aM’(1)6e, (3.28) 


B, ~ 2a(l 


so that 
M(1) = 1 (3.29a) 


M'(1) = -1 (3.29b) 


For the velocity boundary conditions, 
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q; ~ —eU(1 — 67/2)e, + U-be, = 
2UF(1) (1 — 67/2)e, — UF'(1)ée, 3.30 
so that 
F(1) = —<«/2 (3.3la 
F’(1) = -1 (5.51b 


Using the pressure boundary condition 
Pp, = pi + p,U7(1 — 6) (1 — e) (3.32 


along with the boundary conditions in Eqs. (3.29) and 
(3.31), the 6-momentum equation at the shock, 


[u(Ov/Or) + (z/r) (Ov/O0) + (uv/r)|, = [—(1/pr) X 
(0p/00)]; + [(cB,/p) (uB, — vB,)], (3.33 


is evaluated to show that the second derivative of F(R 
at the shock is 


F’'(1) = —(1/e) + 2(1 — €) — 4Q0(1 + €/2) (3.34 


(4) SOLUTION 


Judging from the analysis of the problem, we con- 
clude that we have a system of two nonlinear ordinary 
differential equations for the dependent variables F(R 
and J/(R)—namely, 


F(F’" — 2F"/R — 2F’/R? + 8F/R® )= 
—(2/R?)QeM(— FM" + F’M + 
1FM'/R—4F'M/R) (4.1 


whereas if Q = O, the equation reduces to Eq. 86 of 
Lighthill's paper on hypersonic flow,’ if proper attention 
is paid to notation, and 


—(M”" — 2M/R?) = (2/R?)Rey(—FM'+ F’'M) (4.2 


This system of differential equations is second-order 
with respect to /(R) and third-order with respect to 
F(R). Five boundary conditions are required. 

The five boundary conditions are: 


M(1) = 1 (4.3 
M’(1) = -1 (4.4 
F(1) = —€«/2 (4.5 
ra) = —]J (4.6 
F”’(1) = —(1/e) + 2(11 — €) — 40(1 + €/2) (4.7 


Because of the complexity of the system of differ- 
ential equations, it is not possible to obtain 1/(R) and 
F(R) as tabulated functions, and a computer is neces 
sary. For the computer solution ¢ is taken to be 1/11. 
Magnetic Reynolds Numbers, Re,,, of 1/5 and 1 are 
considered. Values for the magnetic parameter, U, 
range from 1/10 to about 5. 

As can be seen from the analysis, u, the velocity in 
the direction normal to the shock, can be expressed as 


u = 2UF(R) (1 — 6?/2)/R? (4.8 


so that wu is 0 at the value of R for which F(R) 1s 0. 
Thus, a sphere of this radius, to be called R,, is a body 
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shape which is consistent with the spherical shock 
shape and with the flow conditions in the shock layer. 

However, the numerical solution of the computer for 
values of Q greater than a certain critical Q, Q,,, where 
Q., was about 4.2 for both values of Rey, did not 
yield a flow pattern that could be satisfied by a spherical 
body. 

It would appear that considering Q > Q,, in the nu- 
merical solution corresponds to considering Kemp’s 
series solution for his magnetic parameter, S,, greater 
than 4+. For S, > 4, the coefficient of k*/7(=e*’”) in 
the expression for standoff distance, which is propor- 
tional to [1 — (1/4).S,]'/7, becomes imaginary and loses 
its physical significance. In general, it would be wise 
to note that since the problem just supposes a shock 
front, it is possible that a solution to the problem which 
is correct mathematically may not yield a physical 
body. For these reasons the remainder of the report 
is restricted to Q < Q,, where a physical body does exist. 

The following quantities have been computed from 
the solution of the system of differential equations: 
the standoff distance of the shock; the pressure on the 
body; and the gradient of tangential velocity at the 
body. 

(a) The standoff distance 
standoff distance is the distance between the shock- 


The nondimensionalized 


front and the body divided by the body radius—1.e., 
A = (ct — rr) ty = (1 = R,) i, = 
[1 — R(F = 0)]/R(F = 0) (4.9) 


(b) The pressure on the body—To evaluate the pres- 
sure on the body we again use the 6-momentum equa- 


tion, 


u(ov/Or) + (v/r) (Ov/08) + (uv/r) = 


—(1/pr) (0p/00) + (c/p)B,(uB, — vB,) (4.10) 


Expressing the velocity and magnetic terms by means 
of the nondimensional variables, the following equation 
is obtained 


(1/p,,U2) (0p/00) = [40M(—FM’ + F'M)-0/R‘] + 


[((2FF’’ — F’?)0/eR?] (4.11) 
At the body, F = 0, so 
(1/p,,U?) (0p/00)) = 400(M?F’/R‘), — 
(0/e) (F’?/R*), (4.12) 
Integration shows that 
(Po — fr.)/(1/2)e,,U? = 0((40M?2F’/R4) — 
(F’?/eR?)], (4.13) 


if Py is the pressure on the body at the axis of sym- 


metry. Fore = 1/11, 
Poo — Po = (1/2)p,,U°e?P (4.14) 
where 
[(11F’”/R?) 4 
2 = 4.1: 
I k (40F’M2/R°) |, (4.15) 


(ec) Lhe velocity gradient at the body—The velocity 
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TABLE 1 
Rew = 0 
QO A r V, 
0 0.067 2.65 0.524 
TABLE 2 
Rey = 0.200 
Q A P V 
0.50 0.076 2.65 0.406 
1.00 0.088 2.74 0.320 
2.00 0.131 2.94 0.198 
3.00 0.233 3.42 0.120 
3. 20 0.259 3.38 0.103 
3.60 0.374 3.67 0.074 
4.00 0.852 4.32 0.024 
gradient at the body can be expressed as 
[(1/r) (0v/08)], = (U/c)V (4.16 
where 
V, = —(F’/R?), (4.17 
Fig. 2 compares the flow fields for Q = 0 and for 


Q = 0.5, Rey = 1. Tables 1, 2, and 3 show the be- 
havior of A, P, and V, for varying Rey and Q. Figs. 


3, 4, and 5 also show the behavior of these functions. 


(5) CONCLUSIONS 


The first conclusion that can be drawn from the re- 
sults is that the standoff distance of the shock from the 
body increases with increasing Q. Furthermore, the 
standoff distance is relatively independent of Reéy,. 

Secondly, increasing the magnetic field—i.e., increas- 
ing OQ 
given angle. 
Rey has very little effect on the pressure relief. 


increases the pressure relief on the body at any 
Again, for the range of Rey, investigated, 


The third conclusion that can be drawn is that V, de- 
creases as the value of Q increases but is not very de- 
pendent on Rey. Since Meyer® points out that the 
primary mechanism which serves to reduce the heat 
transfer is the reduction of the gradient of the velocity 
at the outer edge of the boundary layer, the reduction 
of V, with increasing Q gives us a measure of how much 
the heat transfer is decreased when the strength of the 
field is increased. 


(Continued on page 728) 


TABLE 3 
Rey = 1.00 
Q A Pr v, 

0.50 0.076 2.64 0.406 
1.00 0.088 2.72 0.320 
2.00 0.131 2.91 0.198 
3.00 0.233 3.36 0.120 
3.20 0.255 3.24 0.102 
3.40 0.299 3.34 0.088 
3.60 0.366 3.56 0.073 
3.80 0.477 3.96 0.057 
4.00 0.754 4.77 


77 0.034 
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Stagnation Temperature Control in 


Blowdown Wind Tunnels 


J. S. MURPHY,* DUNCAN RANNIE,** ann G. W. TIMSON*** 


Douglas Aircraft Company, Inc., and California Institute of Technology 


SUMMARY 


A design method is presented for limiting the variation of 
stagnation temperature in a blowdown wind tunnel so that the 
resultant Reynolds Number variation of the test stream is held 
within acceptable limits. Temperature control is accomplished 
by an installation of thermal mass material consisting of thin 
wall metal tubes or a system of flat and corrugated metal sheets 
with uniform air passages in the downstream end of the storage 
reservoir 

The method presented has been used to design a thermal mass 
installation for a 1 by 1 ft. supersonic blowdown tunnel. Agree- 
ment between calculated and measured air temperature vari- 


itions is shown to be adequate for design purposes 


SYMBOLS 


2 


A = open cross-sectional area of thermal mass tube, ft.? 

A! = surface area of reservoir, ft.? 

a = cross-sectional area of thermal mass tube material, ft.? 

Cp = specific heat at constant pressure, ft.2/sec.? °R. 

( = specific heat at constant volume, ft.?/sec.? °R. 

( = skin-friction coefficient 

¢, = nondimensional heat-transfer coefficient (Stanton Num- 
ber) for flow through a thermal mass tube = 
H puc, 

D = diameter of reservoir tank, ft. 

E = internal energy, ft.-lbs. 


g = acceleration due to gravity, ft./sec.? 
Grashof Number p*gD% 7, — T,)/u?T a 


G, = 

H = heat-transfer coefficient, ft.-Ibs. per sec. °R. ft.? 

h = height of corrugation of thermal mass material, ft. 

k = coefficient of thermal conductivity, ft.-Ibs./sec. °R. ft 

! = length of thermal mass, ft 

M = massof air, slugs 

M = mass flow through tunnel, slugs/sec. 

m = mass flow through unit area of thermal mass, slugs/sec. 
ft.? 

Nu = Nusselt Number 

P = perimeter of thermal mass tube, ft. 

pb = pressure, lbs. /ft.? 

Q = heat, ft.-lbs 

q = rate of heat transfer, ft.-lbs./sec.ft.? 

R= gasconstant, ft.2/sec.? °R. 

Re = Reynolds Number = 4mA /uP 

T = absolute temperature, °R. 

t = time, sec. 

u = velocity, ft./sec. 

V = volume of reservoir, ft.* 

W = work, ft.-lbs. 

x = distance from inlet of thermal mass, ft. 

Z = weight, lbs 

a = heat-transfer parameter = c,P1/A 

8 = heat capacity parameter = amc,tA/p,a,lCps 

\ = variable of integration 
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Mu = coefficient of viscosity, ib.-sec./ft.* 
é = nondimensional length, x// 

p = density, slugs/ft 

o = Prandtl Number = c,u/k 


= nondimensional time = ¢/tpy or t/tgp 


Subscripts 
air in reservoir 


a = 
atm = atmosphere 

e = air at exit of thermai mass 

f = final 

1 = initial, also entrance of thermal mass 
m = reservoir wall 

o == stilling chamber air 

PU = pump-up 

BD = blowdown 

S = thermal mass material 

l = reference condition 


INTRODUCTION 


NE OF THE PROBLEMS peculiar to blowdown wind 
O tunnels is that of maintaining sufficiently con- 
stant stagnation conditions to allow acquisition of test 
data with minimum corrections for Reynolds Number 
effects. Both temperature and pressure must be main- 
tained constant to ensure a constant Reynolds Number 
and test section static pressure. A uniform stagnation 
pressure is ordinarily maintained by means of a servo- 
controlled throttling valve which compensates for the 
timewise pressure variation in the storage reservoir 
upstream of the valve. 

In order to compensate for the decreasing temper- 
ature resulting from the expansion process in the reser- 
voir, heat must be added to the air. An obvious tech- 
nique for accomplishing this is to place in the tunnel 
circuit a properly designed heater which adds the heat 
required to maintain temperature. This 
approach, however, is rather expensive and, as a result, 
several more economical techniques which limit, rather 
than eliminate, the temperature variation due to the 
reservoir expansion process have been suggested. A 
temperature is ad- 


constant 


limited variation of stagnation 
missible since the influence of changing Reynolds 
Number on tunnel Mach Number is a second-order 
effect and small Reynolds Number changes usually 
produce effects on the flow about models which fall 
within the accuracy of experimental measurements. 
Daniels! advocates the use of thin-walled metal cans 
uniformly distributed but randomly oriented through- 
out the high pressure reservoir as a means of limiting 
the temperature drop during the expansion process. 
During pump-up of the reservoir, heat from the in- 
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coming air is stored in the metal which essentially 
reaches the same temperature as the air at the end of 
pump-up. As the air in the reservoir expands during 
a blow, the metal gives up heat and thereby limits the 
temperature drop to an acceptable value. With this 
technique, a temperature drop of 30°F. has been meas- 
ured! during a 35-sec. blowdown from a 9,200-ft.* 
reservoir filled with 157,310 Ibs. of tin cans with air 
initially at 145 psia and 70°F. Disadvantages of this 
method are that it is difficult to design the thermal mass 
installation for arbitrary volume and pressure without 
some recourse to experiment, due to lack of knowledge 
of the magnitude of the heat-transfer coefficient be- 
tween the thermal material and air and the large pres- 
sure drop across the installation. 

Another method of limiting stagnation temperature 
drop is to install thermal mass material consisting of 
metal tubes or a system of flat and corrugated metal 
sheets with uniform passages for the air at one end of 
the reservoir.” Such a thermal mass functions in 
somewhat the same manner as the uniformly dis- 
tributed mass—1.e., it is heated during fill-up and, in 
turn, heats the air passing over it during a_ blow. 
Advantages of this type of installation are that it is 
amenable to calculation based on known thermody- 
namic results and, in addition, pressure drop across the 
thermal mass is negligibly small. It is the purpose of 
this paper to present a complete analysis of a reservoir- 
thermal mass installation of this type and compare 
the calculated results with experiment. 


EQUATIONS GOVERNING THE PumMp-UP AND BLOWDOWN 


PROCESSES 


Consider the system shown in Fig. |. During pump- 
up, air from the compressor flows over the thermal 
mass and enters the storage reservoir. A back pres- 
sure valve is located downstream of the compressor 
to insure constant mass flow during pump-up. It is 
assumed that all entering air passes through the 
thermal mass and is stored in the reservoir. During 
blowdown, air leaves the storage reservoir, passes 
through the thermal mass and the throttling valve into 
the stilling chamber of the tunnel. Neglecting the 
effect of temperature variation of the air leaving the 
thermal mass, it follows that the throttling valve en- 
sures constant mass flow during blowdown. 
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It is assumed in the analysis that the mass flow is q 
constant, positive quantity during pump-up and 4 
constant negative quantity during blowdown; that 
heat transfer takes place by forced convection betweer 
the air and the thermal mass and by free convection 
between the air and the walls of the storage reservoir 
during both the pump-up and blowdown processes 
In addition, it is assumed that heat transfer from the 
reservoir walls to the atmosphere occurs by free con 
vection only. In order to analyze the performance 
of the system, it is necessary to consider a number of 
complete cycles. One cycle is defined as consisting oj 
a pump-up followed by a pause time of variable length 
after the required pressure is reached in the reservoir 
and then a blowdown or tunnel run. The state of 
the system after a cycle is not necessarily the same as 
the state at the end of the previous cycle, since all 
energy put into the system during pump-up may not be 
expended during blowdown. In the analysis presented 
below, the thermal mass and air storage portions of 


the system are treated separately. 


Pump-Up 


Consideration of a single element of the thermal 
mass consisting of an air passage surrounded by thin- 
walled metal leads, by means of a heat balance, to the 
following equations describing the heat-transfer process 
between the air and the thermal mass :° 


Amc,(OT/dx) = —qP ( 


PsCpAs(OT s/t) = GP + kgas(0?T s/0x") (2 


g = ¢,mc,(T — Ts) (3 


It is assumed that m = pu is constant through the 
regenerator. This will be approximately correct if 
the volume of air in the regenerator is small compared 
with the reservoir volume. For a thermal mass dis- 
tributed through the entire reservoir, m varies with 
position and is zero at the far end of the reservoir. 

It may be shown that the conduction term on the 
right-hand side of Eq. (2) may be dropped because it is 
negligible compared to the other terms in the equation. 
Introducing a dimensionless distance from the inlet, 
& = x//, and a dimensionless time, 7 = ¢/tpy, and de- 
fining the dimensionless parameters 
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Schematic of blowdown wind tunnel. 
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" £°” gota 2am 4) 
Eqs. (1)-(3) may be written in the form 
(OT /0&) + a(T — Ts) = 0 (5) 
(OT ;/O07r) — B(T — Ts) = 0 (6) 
The boundary conditions are 
T(0, 7) Ti(r); Ts(&, 0) = T5,(é) (7) 


This system of equations has been solved by Rannie, 
using the Laplace transform method. The solution is 


T= a Ts(& — Ne P~* [py (20V aBrdA)dd + 
J0 


V/ abt /r I, (2 V aBer)dX + « 


V2ax | , 

Verification that Eqs. (8) and (9) are the solutions 
of Eqs. (5) and (6) with the boundary conditions (7) 
may be made by direct substitution using relations (10). 
Note that this solution is not valid for the case of vari- 
able mass flow through the thermal mass since a and 8 
are assumed to be constants. However, the solution 
is valid for a time-varying inlet temperature 7°;(7) and 
any initial distribution of temperature in the thermal 
mass 7°s;(&). 

Air enters the storage reservoir from the thermal 
mass at the temperature 7,,(¢) at the constant rate 
M(lr> 0). Applying the first law of thermodynamics 
to all the air which will be contained in the reservoir 
after pump-up, including that initially in the tank, one 


obtains 


Q=W+AE (13) 


| 


The work done by the gas, IV’, during pump-up is nega- 


tive and of amount 


W=-M [v. p.)dt 


The initial and final internal energies of the air are, 


(14) 


respectively, 


en eee [ rant 
Ey = ae (Mas + Mt) 


(15) 


Hence, making use of the equation of state for a per- 
fect gas, the first law becomes 


t 
Q= -—c,M [ T dt + ¢,T, X 


) 
(Mai + Mt) il CeMaiT at (16) 


Now Q is the amount of heat added to the air during 
« 5 
pump-up and is ordinarily negative since the stored air 





EMPE 





tm?) (5° — 


RATURE CONTROL 


and 


i | Ts(€ — Ae 8" V/ asr/r X 
J0 


Lew ‘aBrr)dy + B } T(r — Ae ™ x 


J 0 


I (2 V/ aBtr)dr +e ”’ Ts (fg) Q) 


The functions /) and /; appearing in these equations 
are modified Bessel functions with the properties‘ 


dJyj/dx = 1,;; dl, dx Ing — (1/2 10 
for small x, 

I(x) > (1/n!)? (x 2)! / 

I(x) > [1/n'(n + 1)!] (wv /2 \ 


and for large x, an asymptotic expansion is* 


tm? | . 


im?) ... }(2n — 1) 


n!(Sx) 


is usually at a higher temperature than the reservoit 
walls. Assuming that the heat transfer between the 
air and the walls takes place by free convection, the 
heat-transfer coefficient may be approximated by 


H = N,k/D = 0.109(k/D) (G,c)'” (17) 


This formula has been confirmed experimentally for tree 
convection from a vertical plate for a range of Grashof 
Numbers appreciably lower than is applicable to the 
problem at hand. However, it is used here for lack of 
a better approximation. 

Then the rate of heat transfer is 


Ga-m = 0.109(k D) [p,2gD* (T. Tn)o/tte?Te)'" X 
ig, ~ 7,3 ti) 
But 
o = Cot/k; pa = PakT; 1 (M,,+ M 
Bo = wi(T,/7 (19) 
so that Eq. (18) becomes 
Ja—m = 0.109¢,[(M,; + Mt) VP? x 


(guy T\0°)' 47 FY" «(oe 


The total quantity of heat transferred to the air during 


QO = -{ qg at (21) 
0 


Hence, the first law becomes 


pump-up is 


—0.109¢,A} (guy T,0?)'’* X 


t 
J, [(Mar + Mt)/VP? (Ty — Tn)? dt 
0 


t 
— cit f Tdt + ¢,Ta(Mai + Mt) — ¢,Me:Tor (22 
0 
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Differentiating Eq. (22) one obtains (with proper ac- 
count for the sign of the term g,— ») 


~0.109¢pA'(gui/Tyo?)'* (Ma; + Mt)/VPO X 
IT. ae Fat’ g (1, aa la i al j ies | = 
—¢pMT, + ¢TaM + ¢,(Ma: + Mt) (T,/dt) (23) 


This equation neglects heat transfer from any entrance 
piping between the compressor discharge and the stor- 
age reservoir inlet. This loss may be appreciable, as 
will be discussed later. 

In order to determine the temperature of the storage 
reservoir walls, the assumption is made that the tem- 
perature gradient within the walls may be neglected. 
An approximate calculation indicates that the resistance 
to heat flow at the inner air-metal face is at least 100 
times greater than the resistance within the metal 
walls; hence, the assumption of uniform temperatures 
within the wall appears justified. Thus, writing a 
heat balance equation across the wall, one obtains 


ee A Oe dt) _ da ~~ Wm atm (24) 


Ja—m iS given by Eq. (20) above, and if heat transfer 
from the wall to the atmosphere occurs by free con- 
vection only, it is given by the relation 


9 Ae Bd 9\1/3 
dm atm = 0.1096 p( Patm™ Zbatm Tomo") x 
AT, — Tatm)’* (25) 


Hence, Eq. (24) becomes 


dT Ma: + uy . ( Gur ) . 
a _ «a = 0.109c » ae ~ A} 
" dt ( if T 0? x 


A ‘ 2s — ya a a atm asi 
IT. — Tn |*”* 5 ) — 0.109¢, (° m = ) x 


. 7 
Te “a Ta! OL atm 


| . 14/3 fe 


Ail™ T atm) 
A 7 nn 7 atm | | qr rad 
| 7 mm ¥ oti 


j (26) 


Eqs. (23) and (26) may be solved simultaneously to 
give storage air temperature and reservoir wall tem- 
perature as functions of time. The boundary condi- 


tions for these equations are 
T,(0) = Tn: and T7,(0) = Ta; 27) 


The air pressure is given by the relation 
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be = (R/V) (Mai + Mt)T, (28 


Blowdown 


During blowdown, the basic equations defining the 
state of the air in the reservoir are similar to those 
used in pump-up [see Eqs. (23) and (26) above], except 
for substitution of 7, for 7, if it is again assumed that 
heat transfer between the air and reservoir walls occurs 
by free convection. As discussed later, this assumption 


appears questionable in the light of experimental 
results. 

Likewise, the equations for the thermal mass soly- 
tion are the same as Eqs. (8) and (9) except for substi- 
tution of 7; for T;. 
is taken as the result from the previous pump-up, with 


the orgin for & displaced to the other end of the thermal 


The value of 7's;(&) for blowdown 


mass. The state of the system including the air in the 
storage reservoir and the temperature distribution 
in the thermal mass at the end of a blowdown is taken 
as the initial condition for the following cycle. 

In case there is a pause time after pump-up and prior 
to blowdown, heat transfer takes place between the 
stored air and the reservoir walls and between the walls 


and the atmosphere by free convection. The equations 


governing this process are 
1/3 2/3 
J M,\* , 
0.109¢,A! () ( ‘) uF)" se 
T\0? ] 
Fe Fs) i a 
IT. T. = —o¢,M, (29) 


—~ 


dT, M,\*" / gu, \'” 
‘omLn —— = 0.109, ( ‘) (; x 
" dt r\y T.0 


1/3 (Ta — Tn) 
IT. — Tm| 


9 1/3 al qn 

P< tmZ Me tm ore on ¢ ( 7 T, tm) 

( i ) A’ |Tn — Tamm |”? 7——"™_ (30) 
a1 atm 7 nn” 7 atm 


AUT, — Tn — 0.109 c, X 


where /, is the mass of air contained in the reservoir 
during the pause time. 


APPLICATION OF THE ANALYSIS TO A PARTICULAR 
INSTALLATION 


Calculated Results 


Using the above analysis, a thermal mass installation 
was designed for a 1 by 1 ft. blowdown wind tunnel 
which has a reservoir with the following characteristics * 


A’ = 3,586 ft.2. Z,, = 82,300 Ibs. 
V = 8,000 ft. p, = 19,920 Ibs./ft.? 
D = 9.5 it. M,; = 18.65 slugs 


Using cp, = 6,000 ft.?/sec.2 °R., C, = 4,285 ft.?/sec.? 
"R., ¢ = 0.72, m = 0.4 XK 10° hs. vec./ft.2, T, = 
560°R., Tatm = 530°R., warm = 0.374 X 10-* Ibs.sec. 
ft.2, and parm = 0.0023 slugs/ft.*, Eqs. (23) and (26) 
governing the process in the storage reservoir become 
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(28 
ee = (iv. = fw 4. These results were obtained by numerical solution 
~20.7( Mai + Mt)? |T. — Tn |*”° = , me nce 
“aie - ’ - T, — Tm of the equation derived above with an IBM 701 Elec- 
tt dT, ; tronic Calculator. The distribution of temperature in 
g the rT IRF ; 4.985T ‘ “ae : 
ir © | 6,000 WT. + 4,285 Mai + Mt) dt + 4,285T,M (31) the thermal mass material at the beginning and end of 
10Se ‘ ° - ° — . . . 
. : pump-up is shown in Fig. 3. The distribution at the 
XCept r 5 ° ae : 
Pu “79 x 108 dT, ~ 90.7(M,, + In?” IT, —T,,|*2 x start of cycle #1 is that measured initially in the thermal 
77 20.7(M, ] a m 
dt = mass in the experiment to be described below. Pause 
CCurs in _ ap ‘ 
ot; 7 = 7.) 7 =n 14 (40 ~ Ten) , time between pump-up and blowdown was taken to be 
tion ; — 145.5 |T,, - : (32 , ‘ 
. a 7 145.5 | Tm — Tatm a sai zero in these calculations. 
ental iF l, . ” vit : ; rages ; Sore” 
: : The timewise variations of reservoir air temperature 
, and the reservoir pressure is related to the temperature and air temperature downstream of the thermal mass 
SOly- hy the e ati ° -" 7 
hat by the equation are shown in Fig. 4. Note that the calculated results 
ISU- . , , ; : : 20% 
ne Po = 0.214(M,; + Mt)T, (33) indicate a temperature drop of approximately 27°F. 
ia : : for a 40-sec. run at 1J = 1.0. Another interesting 
with The installed thermal mass consists of alternate layers Enntecn of ther cthriieand cnuslee te that Gx ences dete 
mal of flat and corrugated sheets, 5 ft. in length with a total not reach an equilibrium state until several cycles have 
the installed weight of 9,850 Ibs. The total cost of the in- been completed 
tion stalled thermal mass was $2,364, or 24 cents per Ib., ‘ 
ken which was broken down into a material cost of $1,354 Experimental Results 
and a labor cost of $1,010. A photograph of the in- es , 
rior war wae i ee Measured values of air temperature in the reservoir 
stallation is shown in Fig. 2. The corrugations are : 
the : rs . ; . ; and at the exit of the thermal mass during a 40-sec. 
in the form of sine waves with a height to width ratio ge a 
alls ‘ : cgeen nae ue es M = 1.0 blowdown are shown in Fig. 5. These meas 
—— of 1/4 and a height of 5/Sin. Thus, the ratio of open i note ' : 
ons : : : - urements were made with 36-gage Iron-Constantan 
° area to perimeter for an air passage is is: 
thermocouples. Measured values of air temperature 
A/P = h/4.28 (34) and pressure in the reservoir during pump-up are pre- 
; ; sented in Fig. 6. The accuracy of the temperature 
The total open area for flow through the mass as in- measurements is +1°F., whereas the pressure was 
stalled is 67.2 ft.2.. The Reynolds Number of the flow measured with an accuracy of +1.5 psia. 
9) through an air passage 1s Additional measurements were made at test section 
Re = (4mA/pP) = Mach Numbers of 0.2 and 3.5. The only unusual feature 
(4g /67.2u) (A/P) = 1,900.17 (35) shown by these measurements occurs in the case of 
M = 3.5 operation, where the stilling chamber was 
ss aiamenens 39 : ' : eee 
or a temperature of 532°R. filled rapidly to a pressure of about 95 psia in order 
P > € “SS ‘livers 0.028 : é eee 
During pump up the compressor delivers ).02S4 to start the tunnel. Due to the rapid stilling chamber 
slug ‘sec. and during blowdown, for the experiment to fill-up process, a thermal pulse was measured, as 
be described later, V 1.943 slugs/sec. ae. shown in Fig. 7. The occurrence of this thermal 
. i a i . 4 y r _— 
) during pump-up, Re = 54, and during bl »wdown, Re % transient is indicative of the fact that the fill-up process 
3,700, so that the flow through the passages is laminar is not isothermal. It is expected that this effect would 
r | for pump-up and turbulent for blowdown. Using the 
| relation 
120 
Cr. = 16/Re (36) ———— START OF PUMP-UP | 
———=— END OF PUMP-UP 
the friction coefficient during pump-up 1s 
1OOF - 
‘sii wit CYCLES 
Crpy = 0.296 (37) No | 
NO. 2 — 
and applying the Reynolds analogy, the heat-transfer 80 er | 
coefficient becomes aoa e | 
~ ) q ~ ° 
| Cy & c,/2 = 0.148 (38) F — 
" CYCLES | 
| 








' For blowdown, experimental results’ for turbulent flow No. | 
| in pipes indicate that : ge | 
Crep = 0.0102, cy ~ c,/2 = 0.0051 — (39) ” | 

| Hence, from Eq. (4), 
20} | 
apy = 60.8, agp = 2.1, Bry = 61, 
Bap = 1.092 (40) 

| fora 5,400-sec. pump-up period and a 40-sec. blowdown. 7 fe) 2 3 4 ) 
Calculated results of the performance of the system DISTANCE ALONG THERMAL MASS-FT 


during three complete cycles are shown in Figs. 3 and Fic. 3. Calculated distribution of temperature in thermal mass 
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Fic. 7. Stagnation temperature variation during a 40-sec. JJ = 


3.5 blowdown. 


be more pronounced at higher Mach Numbers where 
a buildup to higher pressure levels is required for 
starting the tunnel. 


DISCUSSION 


Comparison of the calculated and measured variations 
in air properties in the reservoir and at the exit of the 
thermal mass may be made by referring to Figs. 5 and 
6. Agreement between measured and calculated air 
temperature at the exit is considered to be satisfactory. 
Since determination of this variation is the primary 
goal of the analysis, it follows that the method pre- 
sented may be expected to yield adequate accuracy for 
the design of a thermal mass installation of the type 
considered. However, further examination of the re- 
sults shows that certain portions of the analysis do not 
agree well with experiment, so that improvement of the 
analysis is possible. 

There are two discrepancies between calculated and 
experimental results which will be discussed. These 
are (a) the variation of reservoir air temperature 
during pump-up, in particular during the pump-up 
portion of the first cycle, and (b) the variation of reser- 
voir temperature during the blowdown. As shown in 
Fig. 6(a), the calculated reservoir air temperature 
reaches a peak after about 400 sec. of pumping and then 
decreases to a value about 10°F. greater than the com- 
pressor discharge temperature, which was S7°F.  Ex- 
periment does not show the occurrence of the peak 
temperature and indicates a final reservoir temperature 
about equal to the compressor discharge temperature. 
In the experimental setup, a 30-ft. long, 6-in. diameter 
steel pipe, weighing 19 Ibs./ft., connected the com- 
pressor discharge to the storage reservoir. A rough 
calculation indicates that a heat loss of about 5,300 
3.t.u./hr. occurs to the atmosphere from this entrance 
pipe during pump-up. Since the heat capacity of the 
air flowing through the pipe is about 790 B.t.u./hr. °F., 
it follows that, on the average, the air entering the stor- 
age reservoir has a temperature which is about 6.7°F. 
less than the compressor discharge temperature due 
to this loss. In addition, about 1,980 B.t.u. is re- 
quired to heat this pipe from its initial temperature 
of 70°F. at the beginning of the first pump-up to an 


equilibrium temperature. This results in an addi 
tional drop of the air temperature entering the storage 
reservoir of about 2.5°F. Hence, the air entering the 
reservoir is about 9°F. cooler than the compressor 
discharge temperature. Therefore, it is possible that 
better agreement between experiment and the calcula- 
tions would result by including the heat loss from the 
entry pipe and the heating of this pipe in the calcula- 
tions. The latter effect probably would decrease the 
magnitude of the temperature peak shown in Fig. 6(a). 

The second discrepancy noted above—i.e., the vari 
ation of reservoir air temperature during blowdown 
is such that the calculated rate of heat addition to the 
air from the reservoir walls is less than the measured 
rate. As noted previously, it was assumed in the calcu 
lations that this heat transfer occurs by free convection. 
Consideration of the rate of air withdrawal from the 
reservoir during blowdown indicates that the flow ve 
locity is sufficiently great so that this heat transfer 
actually occurs by forced convection at a rate about 
eight times greater than the free convection rate used 
in the calculations. If this correction were included 
in the analysis, the calculated air temperature gradient 
in the reservoir would agree better with experiment. 
It is evident from a study of Fig. 5 (cycle 1) that both 
of these improvements would yield improved accuracy 
in the predicted air temperature variation at the ther 
mal mass exit. However, it is felt that the accuracy 
of the existing analysis is probably sufficient for design 
purposes. 

The drop in stagnation temperature found for the 
above system results in an increase in Reynolds Number 
with time as shown in Fig. 8, since stagnation pressure 
is controlled at a constant value during a blowdown. 
The effect of this Reynolds Number variation on the 
data measured on models in the tunnel may be found 
by a consideration of its effect on skin friction, transi 
tion, and separation. 

In laminar flow, skin friction is inversely proportional 
to the square root of the Reynolds Number. Hence, 
for a model with completely laminar flow, the skin fric- 
tion at the end of a 40-sec. run would be 3.8 per cent 


60 
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Fic. 8. Timewise variation of Reynolds Number per foot during 
a 40-sec. 1J = 1.0 blowdown 
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less than at the start of the run. If the run were 
only 20 sec. long, the decrease in skin friction would be 
1.9 per cent. It is believed that drag measurements 
may be made to about 1/2 per cent accuracy, so that, 
for completely laminar flow over the model, a correction 
for Reynolds Number variation during a run probably 
is necessary. 

In turbulent flow, skin friction is inversely propor- 
tional to the 1/5 power of the Reynolds Number. 
Hence, in a 40-sec. J = 1.0 run, the decrease in skin 
friction of a model with completely turbulent flow is 
only 1.8 per cent. 

The effect of the variation in Reynolds Number on 
transition may be shown by a consideration of measured 
values of transition Reynolds Number on a cone as de- 
termined in several wind tunnels, which indicate that 
transition occurs at about Re = 6.7 X 10° at WM = 1.0. 
The resultant shift in the x-location of transition due 
to the above Reynolds Number variation with time is 
shown in Fig. 9. Transition moves forward about 1.1 
in. during a 40-sec. run. Experimental error in the 
measurement of transition location is approximately 
0.1 in. to 0.2 in. so that the forward movement of transi- 
tion should be easily measureable in a maximum 
length run. 

The strength of a shock wave necessary to induce 
separation of a boundary layer varies as Re” where 
n = —(1/4) at most for laminar boundary layers* and 
n = —(1/5) to 0 for turbulent boundary layers.’ If 
the Reynolds Number variation with run time shown in 
Fig. 9 is taken to be typical for Mach Numbers greater 
then unity, then the effect of this variation on shock 
strength required to induce separation is either negli- 
gible, or, in a marginal case, a shock wave which is 
about 2 per cent weaker than that required to induce 
separation at the beginning of a run might induce sep- 
aration at the end of a 40-sec. run. 

Since all of the above influences of Reynolds Number 
change are related primarily to force runs, it is apparent 
that the ability to take data for a polar within 10 to 20 
sec. will minimize any corrections for variable Rey- 


nolds Number. Reynolds Number influence on pres- 
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Timewise variation of transition location on a cone 
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sure distribution is a second-order effect, so that no 
correction for Re change or pressure measurements 
should be necessary. 

Additional calculations related to a 
installation for a 6 by 6 ft. blowdown tunnel have 


thermal mass 
shown that an installation similar to that discussed 
above may be designed which gives adequate perform- 
ance for the Mach Number range 0.2 to 5 (see reference 
10). A stagnation temperature of 660°R. may be ob- 
tained for Mach 5.0 operation (to avoid air condensa- 
tion) by regulating the air temperature at the com- 
pressor exhaust to approximately 720°R. Pause time 
after pump-up should be avoided in order to eliminate 
heat loss from the reservoir to the atmosphere and re- 
sultant increases in temperature drop during a run, 
especially at Mach 5.0. Operation under low at 
mospheric temperatures is slightly more critical with 
regard to temperature drop during a run than operation 
at normal atmospheric temperature (530°R.). Oper- 
ation at near maximum dynamic pressure at any 
Mach Number is not critical with regard to thermal 


mass performance. 
CONCLUSIONS 

An analysis of a complete storage reservoir-thermal 
mass installation for a blowdown wind tunnel on the 
basis of known thermodynamic results has been shown 
to agree adequately with experimentally measured air 
temperature variations at the exit of the thermal mass. 
Hence, the method is suitable for design purpose, even 
though some discrepancies exist between the calcula- 
tions and experiment in the reservoir portion of the 
system. 
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Vibrations of Thin Cylindrical Shells Analyzed 
by Means of Donnell-Type Equations’ 


YI-YUAN YU* 
Polytechnic Institute of Brooklyn 


SUMMARY 


The Donnell-type equations of thin cylindrical shells are char 
acterized by being in a completely uncoupled form such that one 
equation of the set involves the radial displacement component 
was the only dependent variable and each of the remaining 
equations relates one of the other displacement components to 
w. Derivation is first given for Donnell-type dynamic equations 
including and excluding, respectively, the effects of transverse 
shear deformation and rotational inertia. The equations are then 
applied in a convenient and straightforward manner to the analy- 
sis of free vibrations of infinitely long cylindrical shells, with 
comprehensive results given for the frequency, modes, and phase 
and group velocities. Results on the frequency and modes are 
also immediately applicable to finite shells with simply supported 
edges. The cases of finite shells with clamped and _ flexibly 
supported edges are finally discussed by further making use of 
method in conjunction with the Donnell-type 


the Galerkin 


equations. 


SYMBOLS 


X,S,2 = coordinates in longitudinal, circumfer- 
ential, and radially inward directions 


as shown in Fig. 1 


r = coordinate in radially outward direction, 
related to z by the relation r = a — 2 

Cx, Ce, Cr = normal stresses 

a a = shearing stresses 

“a = normal strains 

Yr» Yzr> ¥ = shearing strains 

i, 0, W = displacements in x, s, 2 directions 

M,, M, = bending moments 

Ms, Mex = twisting moments 

N,, N., Nzs, Nzz = membrane forces 

0,, 0. = transverse shearing forces 

u,v, W = displacements of middle surface of shell 


in x, Ss, directions 
b, ¢ = changes in slope of normal to middle sur- 
face of shell in x, s directions 
amplitudes of u, v, w, ¥, 


time 


t = 
E = Young’s modulus 

m = Poisson’s ratio 

G = shear modulus 

D = flexural rigidity of shell = Eh*/12(1 — yp?) 
p = density of shell material 

k = shear constant 

K = «(1 — p)/2 

l = length of shell 

a = radius of shell 

h = thickness of shell 


Presented at the Missile Structures Session, Twenty-Sixth 
Annual Meeting, IAS, New York, January 27-30, 1958. 

+ This research was supported by the Air Force Office of Scien- 
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* Professor, Department of Mechanical Engineering. 
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k = h?/12a? 

L = wavelength 

m = number of circumferential waves 

n = number of longitudinal half-waves 

r = 2ra/L for infinite shells, or nza for 
finite shells 

Y, = m? + (nza//l)? 

w = circular frequency 

Q2 = (1 — pe E pa*w* 

( = phase velocity 

( = E p\ l Ti ats 

r a aia 

oy = group velocity 

F,G = functions defined in Eq. (5.9) or Eq 
(6.5 

VW = function defined by Eq. (9.4 

v2 = (Oo? Ox?) 4 (Oo? Os* 

0?/oT? = [(1 — u?)/E) pa*(d?/or? 

SR = abbreviation for effects of transverse 
shear deformation and rotational in 
ertia 

1 subscript indicating effect of rotational in- 
ertia 


(1) INTRODUCTION 


Sx the thin cylindrical shell is one of the basic 
elements commonly used in airframe and other 
structures, a knowledge of its vibration characteristics 
has been of practical importance as well as of theoreti- 
cal interest. The importance of the vibration problem 
of thin cylindrical shells has further been recognized 
recently in connection with the design of engine cases 
of jet engines, as it was found that vibration of the en- 
gine case resulting from excitations from primary burn- 
ers and afterburners could be so severe as to cause fail- 
ure of the case. A similar situation could arise in a 
rocket engine, which often has an engine case in the 
form of a cylindrical shell. The problem of vibration 
of the cylindrical shell is therefore becoming increas- 
ingly important, and it is our purpose here to present 
an account of the latest development and results on 
this subject, including those based on the theory which 
takes into consideration the effects of transverse shear 
deformation and rotational inertia. (These effects will 
hereafter be referred to as SR.) It has been known 
for a long time in the case of beams that these effects 
play an important role in vibrations involving high 
frequencies. 

The method of approach used here in obtaining the 
results is associated with the so-called Donnell-type 
equations of cylindrical shells, which are of such a 
completely uncoupled form that one equation of the set 
involves the radial displacement component w as the 
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Fic. 1. 


only dependent variable and each of the remaining 
equations relates one of the other displacement com- 
ponents tow. This type of equation was first adopted 
by Donnell! in the investigation of the buckling of 
His original static equations have 
A brief review of their applica- 


cylindrical shells. 
since been widely used. 
tions may be found in reference 2; we mention here 
only N. J. Hoff* and his coworkers, who have made 
extensive use of these equations in the bending as well 
as buckling analysis of cylindrical shells. 


Donnell’s original static equations were derived by 
assuming that all strain expressions except that for the 
circumferential membrane strain are the same as those 
used in the classical thin-plate theories. In the ex- 
pression for the circumferential membrane strain an 
additional term is included to account for the effect of 
the radial displacement of the shell. The Donnell-type 
equations used in this paper are of course dynamic 
equations, and they do not involve Donnell’s simplify- 
ing assumptions for the strain expressions. Further- 
more, the effects of SR are also incorporated in one of 
the two sets of equations used. These effects are 
readily suppressed to obtain the other set of equations, 
which thus constitutes the exact counterpart of the 
first set. A derivation of these equations is given in 
some detail in Section (II). 


The two sets of Donnell-type equations, which include 
and exclude SR respectively, are then applied to the 
problem of vibration of infinitely long cylindrical shells 
in Section (III). In solving this problem, we note the 
advantage of employing Donnell-type equations in that 
they readily yield compact expressions for the fre- 
quency, modes, and phase and group velocities. De- 
tailed numerical results have also been obtained, and 
some of these are presented here in graphical form. 
More features of both the theory including SR and the 
theory excluding SR are revealed here than those given 
in previous publications, and the shortcomings of the 
latter theory are clearly indicated. Only nonsymmetri- 
cal vibrations are discussed in this paper; axially sym- 
metrical motions of cylindrical shells have been dis- 
cussed by Herrmann and Mirsky,‘ who established 
the importance of SR in the shell theory by showing 
the excellent agreement between results obtained for 
the phase velocity from both the shell theory and the 
three-dimensional elasticity theory in the axially sym- 
metrical case. The problem of propagation of nonsym- 
metrical waves in infinite shells was also recently in- 
vestigated by Mirsky and Herrmann’ on the basis of 
the theory including SR. However, Donnell-type 
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SCIENCES 
equations were not used by these authors. The phag 
velocity was thus computed directly from an equation 
in the form of an unexpanded determinant of the fifth 
order equated to zero, and no other quantities besides 
the phase velocity were considered in their work, pre- 
sumably due to the complexity in computation which 
would have been involved without the use of Donnell 


type equations. 


In Section (IV) a discussion is given of the vibration 
of finite shells. It is shown first that the results on the 
frequency and modes obtained in Section (III) for in- 
finite shells are immediately applicable also to finite 
shells with simply supported edges, when either the 
theory including SR or the theory excluding SR is 
used. This in the latter case has been well known. 
Finite shells with clamped and flexibly supported edges 
are finally discussed on the basis of equations not in- 
cluding SR. For these cases Donnell-type equations 
are used in conjunction with Galerkin’s method, thereby 
illustrating another advantage of these equations. 


In this paper we consider only thin cylindrical shells 
namely, shells whose thickness-radius ratio h/a is so 
small that the parameter k = h?/12a* may be neglected 
in comparison with unity. This requirement is cer- 
tainly fulfilled, for instance, for ordinary cylindrical 
shells, the thickness-radius ratio of which is smaller 
than 1/10. 


(II) DYNAMIC EQUATIONS OF CYLINDRICAL 
SHELLS 


The dynamic equations of cylindrical shells including 
SR as given in reference 6 were derived in the manner 
of Mindlin’ and Herrmann and Mirsky.* | While the 
former author* introduced equations including SR for 
plates, the latter authors further used the same method 
to derive equations for the axially symmetrical motion 
of cylindrical shells. In the following the derivation 
which led to the equations given in reference 6 is 
described. The equations are finally arranged to appear 
in the Donnell type. 


(1) SHELL-STRESS COMPONENTS 


The coordinate system (x, s, z) to which the shell is 
referred is shown in Fig. 1, with z related to 7 in the 
usual cylindrical coordinate system by the expression 
ry = a — 8, where a is the radius of the shell. Since we 
are concerned with the free vibration, the cylindrical 
surfaces of the shell are assumed to be free from ex- 
ternalloading. Thus, 
= %, = 7, =O, = = h/2) (1.1) 


The moments, membrane forces, and_ transverse 
shearing forces, all per unit length, are defined in the 
usual manner: 


* According to Mindlin,’? the method was first used by J. 
Boussinesq to derive plate equations. 
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which will be referred to as the shell-stress components. 


(2) RELATION BETWEEN SHELL-STRESSES AND SHELL- 
DISPLACEMENTS 


We make use of Hooke’s law in the following form: 


oe = [E/() — p*)) tee + es) + 
IM/(1 — pw) Jo 

= (£/(1 — g°)| te + wee) + (2.1 
lu/(1 — pw) Io 

te = Gren te = Gite ta = Gy 


1O O O 
a D( Ws ee *’, M 
aox Ox Os 
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where the strains are given by 


€, = OUW/OX, €, = (a/r) (O0/Os (@/r 

Yrs = (a/r) (OU/Os) + (O0/Ox), (99 
Yer = (O0U/Or) — (O@/Ox aii 
_:; = (O7/Or) — (@/r) — (a/r) (OW /Os 


The displacements are next assumed in the form 


“= u(x, s,t) — ep(x, s,t 
? = v(x, Sf 7O(x, S, lt 2.3 
a w(x, Ss, f 


where #, v7, w are the displacements of the middle. sur 


face of the shell in the x, s, z directions, and Y and @ 
are the changes of slope of the normal to the middle 
surface in the x, s directions, respectively. These five 
quantities are used to describe the motion of the shell 
and may therefore be called the shell-displacement 
components. 

The shell-stress components can now be related to 
the shell-displacement components through the suc 
to (2.3) in (1.2); the 


cessive substitution of Eqs. (2.1 
results after integration are 


1 Ov WwW oy Od 
=» Di. Pr 4 
a Os a Ox Os 


| Ou oy ) 


l 1 Ov Oo Op ] 
M., = Da u) ( + v + ), M,, = —-D(1 - u) (- : 
2 a Ox Os Ox 2 aos Os On 
‘ kh Ou Ov w h? ow ; kh Ou O% wW h? Om 
Nz = = +a a 4 , Ny = -ta— + --- 2.4 
1 — p* \Ox Os a 12a Ox 1 — pe Ox Os a 12a Os 
— Eh (~ ov h? oe) Eh (= Oz h? *) 
7 "211 + w) \Os © Ox 12a Ox/’ 21+) \ds ox 12a 9s 


‘ Eh ( ow ’) 
ee Bt ox 


In carrying out the integrations in Eqs. (1.2), the fol- 
lowing rules have been used. (a) Integrals involving 
o, are neglected. (b) Results for the transverse shear- 
ing forces Q, and Q, are multiplied by a shear constant 
x. (c) The factors a/r appearing in Eqs. (2.2) are ex- 
panded in series form as 


(z° 


a/r = a/(a — 2) = 1 + (2/a) + (2°/a*) + 


which is then introduced in the integrals in Eqs. (1.2). 
In the results of integration use is finally made of the 
thin-shell requirement that k = h*/12a°’ < 1. 

Rules (a) and (b) were first used by Mindlin’ in the 
derivation of the corresponding relations in his plate 
theory, and then by Herrmann and Mirsky* in the 
study of the axially symmetrical motion of cylindrical 
shells. The latter authors further used a third rule 
that logarithmic terms are expanded in series of as- 
cending powers of /:/a and terms of a degree higher than 
the third are neglected. Although adoption of their 
third rule would lead to the same results in the non- 
symmetrical case, the above rule (c) is clearly less 
arbitrary and therefore preferred. 

The shear constant «x may be determined by either 
of two methods, both of which were first used in the 
plate theory by Mindlin and extended to the shell 


Eh {f ow =? ) 
- - T QD 
2(1 + Lu ( Os a 


theory by Herrmann and Mirsky. In the first method, 
x is chosen in such a manner that for infinitely short 
waves the shell equations yield the frequency of the 
Rayleigh surface waves. This yields the result 


») 


— [(1 — 2u)/2(1 — pw) kj [1 — («/2 


{ 
(1 — x) j}1 


The desired value of x is given by the lowest root of this 
equation, which depends on Poisson's ratio. For u 
0.3, we find 

x = 0.36 


In the second method of determining «x, frequencies 
calculated from the shell equations and from the elas- 
ticity equations for the thickness-shear vibration are 
matched, which yields 


kK = 7/12 = 0.824 


Although the first method of determining « was criti- 
cized by Fliigge,” the value of « = 0.86 had already 
been used in the numerical calculations to be discussed 
later, and no change of its value was made, as the dif- 
ference between the results of the two methods is small 
and actually vanishes when uw = 0.176. 
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(3) DyNAMIC EQUATIONS OF CYLINDRICAL SHELLS 


The dynamic equations of cylindrical shells are based 
upon the fundamental equations of motion in the 
theory of elasticity, which, in cylindrical coordinates, 
are 


(O0,/Ox) + (@/r) (Orsz/OsS) + (O7;2/Or) + 
(7,z/7r) p(07*a/Ot?) 
Ox) + (07,;/Or) + (2 
(27,;/r) = p(O070/Odt?) 
(Oa,/Or) + [(o, — o,)/r] + (O7;,/0x) + 
(a/r) (O7s7/OS) — p(0°@/dt?) 


(a/r) (O0;/Os) + (O07; 


Five equations for the shell will be derived from Eqs 
(3.1). The first three are obtained by multiplying the 
latter by v/a, writing a, 3, @ according to Eqs. (2.3), and 
integrating with respect to z over the shell thickness. 
The other two are obtained by multiplying the first 
two of Eqs. (3.1) by sr/a, writing a, 3 according to 
Eqs. (2.3), and integrating. The results are, after mak- 
ing further use of the conditions in Eqs. (1.1), 


CE 
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(ON,/0x) + (ON;,/Os) = ph[(O07u/Ot?) + 
(h2/12a) (0%y/dt2) 


(ON,/Os) + (ON,,/0x) + (Q,/a) = 
ph[(0%/dt2) + (h2/12a) (0°b/dt?) 
(0Q0,/0x) + (0Q,/0s) — (N,/a) = : on 
— ph(d°w/Odt?) — 
(OM,/0x) + (0OM,,/0s) — QO, = 


p(h?/12) [(1/a) (07u/ot?) + (07y/ot?) ] 
(OM,/Os) + (0M,,/0x) — QO; = 
p(h®/12) [(1/a) (0°v/ot?) + (0°¢/odt?)] | 


Eqs. (3.2) were also derived by Mirsky and Herrmann,* 
by means of the energy approach, who further formu- 
lated the appropriate boundary and initial conditions, 
The boundary conditions for a closed shell are such that 
at each of the two edges one member each of the five 
products M/,y, \/,,.¢, N,u, N,v, and Q,;w must be pre- 
scribed. By letting the radius a equal infinity, the last 
three of Eqs. (3.2) are reduced to Mindlin’s results for 
a plate.’ 

Eqs. (3.2) are in terms of the shell-stress components. 
By means of Eqs. (2.4), we may further put these 
equations in terms of the shell-displacement compo- 
nents; thus, 


Ou 1 — p Ou l1+yu Ov pb OW +k (- l—u a4 1— pu? Ou 4b l1—p’ Oy 
_ a — = 
ax? 2 ds? 2 daxds  adx Ox? 2 ds? gE Po cea 
1+ yu 0% 1—pOw Ow 1 Ow (’ — 20° *) 
: -— + ka - — — _ 
2 Oxods 2 Of Os? a Os 2 Ox? Os? 
_fv 1 Ow o 1—p Ov 1— yp OO} 
K _ = k 
( a OS *) E 4 or? 7 E P dt? 
ra) 1 Ov w re) 1 Ov ra) Oo 1—p? Ox 
yu Ou Se + K/( hia :. # SS 
a Ox a Os a? Os a Os Ox Os / f or? 
(3.3) 
0” 1 — »O? oO 1 — oO” l 0" 
(Ss ro +s v Pa v ee 4 
Ox? 2 Os? x? 2 Os? 2 OxOs 
_ fl dw “) 1—p? O% 1—p? dy 
K _ =k k 
(; dx a E Poe) Pape 
1—pow Ov 1d0w l1+4u o%y l—p 0% *) 
k 2 3 ye + 
( . ow ao” Sl Cee Se 
_fv 1 Ow ) —p? Ov 1—p? 0% 
K{— —_— = = — k 
(*. + a Os a EL ’ Or? la ihe Of)? 


where subscripts 1 have been added to the time variable ¢ in all terms affected by rotational inertia so that its effect 


may later be traced. 


When SR are neglected, Eqs. (3.3) may be shown to reduce to 


O7u 1— uO 1+u 0% yp OW 3qy 
— : — — — + ka - — 
Ox? 2 Qs? 2 Oxos a Ox ox? 
l+u Oru 1 — pO” O*v 1 Ow 
2 oOxods 2 ox Os? a Os 
31— uOo%% 
ka ( Pee a 
2 a Ox? 


1— 
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1 0°v 
a Os? 


O7u 
? ot? 
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7 oe) 7 | 
dxds2] 
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y Ou 1 Ov w , (| O*u l1—yp O*u 3—y4 ev 1 Ov 2 O-w os ) 7 
? a Ox* 2a Oxds? © 2a Ox"ds a*® Os a” Os* ” | 


3.4) 


p 


1—p Ow cont 
3.2) E ot- 


-. These differ from the well-known Fliigge equations® in the small & terms only; the latter may be obtained from Eqs. 
3.4) by inserting the additional term [(1 — u)/2a] (0°u/Os*) in the parentheses of the first of these equations, by 
dropping the terms —(1/a) (0%”/0s*) + (1/a*) (Ow/Os) in the parentheses of the second, and by replacing the term 
|/a*) (Ov/Os) by the term w/a‘ in the parentheses of the third. Since k is negligibly small compared with unity for 





ordinary thin cylindrical shells, the difference between Eqs. (3.4) and Fliigge’s equations in the coupled form is 


nn,? 
mu- actually very small. 
ONS, nn - 
hat (4) DONNELL-TyPE DyNaAmiIc EQUATIONS 
five The equations used in our analysis are Donnell-type equations, the first set of which is derived from Eqs. (3.3), by 
aa differentiation and elimination. Details of the derivation are described in reference 2, with the final results of the 
ast equations given by 
for 
rk / oe li-—_ .. k [o> os oO l—-uw.. 
ts |—{—— ~ a°v?) + 1|| —(—— — a’v?) +1 — = a°v?) X 
re LK or;* 2 K oT; ol” 2 
“se 
: d° O*w O° l1—yp , Oo o*w | k { l—u 
a -a°Vv*? + i) —) ( = a* - 1} (1 — gw?*)a? > | —k ; se av. i+tli x 
o7 of oT? 2 Ox Ox K \oT; 2 


3 ‘|| “ | ; . | |  ( - | ‘) | | | 
—— — g?V? — — av? — — @V?) av2w — 2ua + 
or? or? 2 oT: ox? 


O*w O°w O'w l—s 7—3u Ow O'w 
k(1 — p)a®| (38 — 2p) + 2(2 — pu | +4 a‘| + = Q (4.1 
| , Ox'4ds? 6 9) 9 20852 Os} | 


i: ( = -v*] +1 | a | | ~ | iF ( “4 v?) | x 
-— oe = a°“V “ -~ iy ia eS 4 = ow — 
a ye or 2 ° or (LK \o7? 


| 0? l1—yp ,. Oo l—u ] a(S l—yp . oe =) 
m -— pwa* > a? — | + k — ae a- + ywa- — 
oi 2 Ox 4 Os? o1;*\or* 2 Ox Os* 








o 1 — o? 1 - o? ra) O 3 — O° 1 — a) 
| ’ (2+ n)a? = ra |+2 = (5 - a? _— =a ) - 
oT? 2 ox? 2 Os 07,2 \OT? | Ox? 2 Os 
a) 3 — pw . Oo o? | mn 
k 1° +k(1 — p)a'v? >a 1.3 
nat 2 * =) Oy eet” Bs 
'(k/K) [(02/07,2) — a?V2] + 1f y =Ow/On 1.4 
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! ; a =a? ‘) | Il 3 ¥ : | | “ | | 
™ we) 4] ia ay? _ av? | — 
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E Sed 2" a) l—u Rea + 
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When the various terms in Eq. (4.1) are expanded and some of the small k terms neglected according to the rule that 
k is negligibly small compared with unity, we obtain 
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O° ]1 — e . O*w 0” l O°w [ - Oo” 1 — Me ’ 
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Eq. (4.6) is in every respect the same as Eq. (4.1) when the former together with Eqs. (4.2) to (4.5) is applied to the 
analysis of nonsymmetrical vibration of thin cylindrical shells. When the motion is axially symmetrical and thus 
independent of s, Eq. (4.6) will not yield the proper result, but Eq. (4.1) can be reduced to the proper form for the 
special case. Details of the latter reduction may be found in reference 2. 

A second set of Donnell-type equations which do not include SR and are also used in our discussion of the vibration 
problems may similarly be derived from Eqs. (3.4), but it is more convenient to deduce them from Eqs. (4.1) to (4.6 
by putting equal to zero terms involving the factor k/K and the operator 0°/07)°, which represent the influence of 
transverse shear deformation and rotational inertia, respectively. Thus, Eqs. (4.1) to (4.3) become 


( °* l—p v2) Q* 272 4 1) O°-w ( Oo” oe ) (1 ew 4 
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or 2 Ox? Ox* | Os 
Similarly, Eq. (4.4) reduces to Y = Ow/Ox, which is simply the condition that the transverse shear deformation 


vanishes in the x-direction. Eq. (4.5) cannot reduce to such a simple form, because the transverse shear deformation 
in the s-direction involves v in addition to ¢ and w which are present in the equation. When the terms affected by 
SR are eliminated from Eq. (4.6), the following result is obtained: 
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This would be the same as Eq. (4.7) if the small term —2ku(0?/07")a‘V?(0°w/Ox?) in the latter were not now missing. 
Eqs. (4.8) to (4.10) have also been used in our study of nonsy mmetrical vibrations of cylindrical shells. The accu- 
racy of the results is of course not affected by the absence of the negligibly small k terms from Eq. (4.10). 

When the coefficient (7 — 3u)/2 in the last term of Eq. (4.10) is replaced by 2(2 — u), the equation becomes one 
of the Donnell-type Fliigge equations given in reference 6. The other two equations of this last set are again the 
same as Eqs. (4.8) and (4.9). When all & terms are neglected in Eqs. (4.8) and (4.9) and & terms of the fourth and 
sixth order neglected in Eq. (4.10), we further find the equations which constitute the exact dynamic counterpart 
of the original Donnell static equations as were derived more directly in reference 9. A detailed discussion of Donnell 
and Donnell-type equations may be found in reference 2. 
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(II) VIBRATIONS OF INFINITE SHELLS 





(5) VIBRATIONS OF INFINITE SHELLS BASED ON EQUATIONS INCLUDING SR 


xX The shell-displacements are assumed in the following form: 

u = LU cos (Ax/a) cos (ms/a) sin wt 

v = V sin (Ax/a) sin (ms/a) sin wt 

w= W sin (Ax/a) cos (ms/a) sin wl (5. 
(4.6 y = WV cos (Ax/a) cos (ms/a) sin wl 

@ = ® sin (Ax /a) sin (ms/a) sin wt 
O the . ps , . 9 -f . 
where is the number of circumferential waves equal to or greater than one, and A 2ra/l may vary from zero for 
thus a ; 


infinitely long waves to infinity for infinitely short waves. To find the natural frequency and mode of vibration of 
r th ES ow : : ; 7 : = 
. infinitely long shells, these displacements are substituted in Eqs. (4.2) to (4.6). The results are 
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| in which Q = [(1 — u?)/E] pa’w? and Q; denotes the same quantity derived from 0°/0f,;? terms and is otherwise 


entirely the same as 2. Eq. (5.2) is the frequency equation. For any given material and dimensions of the shell 
| and a given value of m, the equation yields five frequencies for any value of \ between zero and infinity. When the 
value of 2 becomes known and is substituted in Eqs. (5.3) to (5.6), the corresponding modes may be calculated 
straightforwardly. It is to be noted that the frequency equation (5.2) is good for any value of m 2 1. A similar 
frequency equation immediately derivable by substituting w from Eq. (5.1) into Eq. (4.1) instead of Eq. (4.6) 1s also 
good for m = 0; namely, for the axially symmetrical case and the torsional case. 
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The problem of vibration is closely associated with that of wave propagation. To find the phase velocity ¢ = 
wa/r, we note that 


Q/d? = [(1 — w*)/E]pa?w?(1/d*) = c?/[E/p(1 — u?)] 


Denoting co? = E/p(1 — p?), T= c2/e9? 
we have Q/r2 = c2/e? = T 
Similarly, we indicate the effect of rotational inertia by writing 2,/A2 = ¢,2/co2 = T;.. In terms of I, Eq. (5.2) may be 


written as 
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(1 —- (1 — p?) 1 1 — 1 /7-3 
# alt -— k 4 ( = 2m? + m')= 0 (5.7 
2 r$ 2 x 2 


As Eq. (5.2) yields five frequencies for any given values of k, m, and X, Eq. (5.7) also has five roots of T in each case. 
In a problem of wave propagation in a dispersive system, another quantity of importance is the group velocity, 
which by definition is the velocity of propagation of a wave packet consisting of waves whose wavelengths do not 
differ much from a certain given value. It may be shown that energy is transmitted by the wave with this instead 
of the phase velocity. The group velocity C is related to the phase velocity c and wavelength L as follows: 
C = c — L(dc/dL) 
In terms of \ this becomes C/co = (c/co) + Ald(c/co)/dd] or C/o = YT + A(dVT/dd) 
Remembering 1/1 = 1/Q/A, we obtain finally = C/co = d+~/Q/dd 5.8) 
It is seen that, if ~/Q is plotted against \, the slope of the resulting curve is equal to C/co. By making use of this 
relation and with +/Q or +/T given implicitly in terms of \ in Eq. (5.2) or (5.7), another expression of the group 
velocity better for the purpose of calculation may further be derived—thus, 
C/eo = (1/VTP) [F(P, d)/G(IP, d)] 5.9) 


where 
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In the functions F and G, k; is now used to trace the effect of rotational inertia and is otherwise entirely the same as &. 
The limiting cases of infinitely long and short waves, for which X = 0 and A = »&, respectively, are of particular 
When A = 0, Eq. (5.2) yields five finite values of 2, and the corresponding values of [ are therefore all 


interest. 
= ], for which case it may be shown from Eqs. (3.3) that the lowest value of 2 is zero and 
; | 


infinity (except when m = 


the lowest value of Tis finite). On the other hand, when A = ~, Eq. (5.7) reduces to 


(T — 1) (PF — [0 — #)/2]; (@/K)P — ke] @ — 1) ie — (Cl — w)/2]j = 0 


which yields immediately 
r= K, (1 — p)/2, (1 — g)/2, 1, 1 (5.10 


and, consequently, Q = Ad, [((1 — w)/2]d2, [C1 — w)/2)A2, A, A? 
The frequencies therefore all increase with \ indefinitely. 
The limiting values of the group velocity in these special cases may also be deduced easily. When A = 0, the frac- 


tion F/G in Eq. (5.9) is always finite. Thus, since all five values of I are infinite in this case, the corresponding values 


of C/cy are all zero. When \ = ©, the values of T are given by Eq. (5.10). Corresponding to the value T = A, 
Eq. (5.9) yields directly (C/co)? = AK. With the other values Tf = 1 andl = (1 — u)/2 substituted in Eq. (5.9), 
The results for (C/co)? 


the fraction F/G becomes of the indeterminate form 0/0, and L’Hospital’s rule must be used. 


turn out to be always the same as for T'; namely, 
(C/ta)* = K, (1 — »)/2, (1 — g)/2, 1, 1 


The group velocity is therefore always finite for infinitely short waves according to equations including SR, as it 


should be. 
(6) VIBRATIONS OF INFINITE SHELLS BASED ON EQUATIONS EXCLUDING SR 


The same procedure which has just been used in solving the vibration problem of infinite shells may still be fol- 
lowed when the analysis is to be based on Eqs. (4.8) to (4.10), which do not include SR, instead of Eqs. (4.2) to (4.6 
The same results however may be obtained more conveniently by putting all terms affected by SR equal to zero in 
Eqs. (5.2) to (5.4), (5.7), and (5.9); namely, all those involving k/AK orthe subscript 1. Thus, these equations reduce to 
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For given values of k, \, and m, the frequency equa- 
tion (6.1) now yields only three instead of five fre- 
quencies. It may be seen that, when \ = 0, the three 
frequencies are finite, the corresponding phase veloci- 
ties are infinite (except when m = 1, for which case 
Eqs. (6.1) and (6.4) yield zero and k(1 — yu)/4, respec- 
tively, for the lowest values of Q and T), and the group 
velocities are zero, just as when equations including 
SR were used. When } = ©, however, the phase 
velocities, according to Eq. (6.4), are given by 


1, (1 — p)/2, 
from which 
Q = d*, [(1 — w)/2]r2, 


The corresponding group velocities, as calculated from 
Eq. (6.5), are given by 
(C/ceo)? = 1, (1 — »)/2, 

which are again all equal to T, as when SR were con- 
sidered. The infinite value of the group velocity for 
infinitely short waves implies that energy would be 
transmitted at an infinite speed by such waves, which 
is physically impossible. This constitutes one of the 
main differences between the results produced, respec- 
tively, by equations including and excluding SR, and 
is, in fact, one of the principal shortcomings of the latter 
equations. 


(7) NUMERICAL RESULTS FOR INFINITE SHELLS 


The results obtained in Eqs. (5.2) to (5.4), (5.7), 
(5.9), and (6.1) to (6.5) for the infinite shell may be ex- 
amined further through numerical calculations, which 
have been carried out for the values h/a = 1/30, 1/10, 
and m = 2,5. There are thus four combinations, and 
in each of them calculations were made for most or all 


of the following values of \:* 


\ = 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500 


which will be seen to cover practically its full range. 
The value of » was taken equal to 0.3 and x equal to 
0.86, corresponding to which A has the value 0.301. 
The quantities calculated include 2, 1/2, WT, C/c, 
U/W, and V/W, which are all dimensionless. 
these, 


Among 
Q, WT, and C/c are proportional to the fre- 
quency, phase velocity, and group velocity, respec- 
tively, and U/W and V’/W represent the mode. Most 
of the calculations were made on an IBM 650 electronic 
digital computer, with the results given in eight sig- 
nificant figures, and the rest of the calculations were 
carried out with the aid of a table of logarithms or a 
slide rule. The numerical results for all the above 
quantities were tabulated in reference 10, but they are 
only partly presented here in graphical form. Thus, 
the graphs of ~/T, plotted versus } for all four cases of 


* Also other values to be mentioned later. 
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IBRATIONS OF 
bk = 1/30, 1/10, and m = 2, 5, are shown in Figs. 2 
to 5, and those of the other quantities as well as fur- 
ther details of »/T for only the case of k = 1/30 and 


m = 2in Figs. 6to 12. 


Phase Velocity 


The graphs of ~/T versus \ in Figs. 2 to 5 represent 
the dispersion curves. The solid and dotted curves 
in these (and all other) Figures indicate the results 
based on equations including SR and those based on 
equations excluding SR, respectively. The dispersion 
curve thus has five branches when SR are considered 
but only three branches when SR are disregarded. The 
limiting values of these branches are clearly the same 
as those deduced previously for the general case. The 
agreement between the solid and dotted curves is seen 


to be excellent up to the value A = 15 in Figs. 2 and 3 
for the two cases of h/a = 1/30 and up to A = 6 in 
Fig. 4 for the case of h/a = 1/10 and m = 2. For 
the case of h/a = 1/10 and m = 5 in Fig. 5 the agree- 


ment is generally good up to \ = 5, but there is always 
some noticeable discrepancy for the lowest branch even 
for \ < 5. Beyond these aforementioned values of \ 
for the various cases, not only the lowest branches of 
the solid and dotted curves show large deviations from 
each other, but also the arrangement of the branches of 
the two types of curves becomes quite different, as 
may be seen from these Figures. Comparison between 
the solid and dotted curves beyond these limits is there- 
fore hardly feasible, except perhaps just for the values 

















of VT in the case of X approaching infinity, which, 
however, was already discussed. 


It is noted that, for both the solid and dotted curves 
in Figs. 2 to 5, the lowest branch in each of the two 
cases m = 2 displays two minimum points, while in the 
two cases m = 5 each has only one. This character- 
istic of the lowest branch in the case of the solid curve 
for which SR are considered was also observed by 
Mirsky and Herrmann.’ Comparing the four Figures 
we also conclude that varying the value of m from 2 to 
5) but keeping /,/a constant mainly affects the left-hand 
half of the graph with the right-hand half remaining 
practically unchanged. On the other hand, varying 
h/a from 1/30 to 1/10 and keeping m constant cause 
only very little change on the left-hand portions of the 
second and third branches of the curves. 


Results for V T based on equations including SR have 
also been reported by Mirsky and Herrmann.’ Good 
agreement was found between their results and the re- 
sults given here, although x was taken equal to 7°/12 
in the former and 0.86 in the latter. Not for the dif- 
ference in the value of x, however, it is believed that 
accurate calculation of WI is much easier by using 
Eq. (5.2) than by using the unexpanded fifth-order 
determinantal equation in reference 5. 


It is seen that in each of Figs. 2 to 5 there are three 
locations where, as a result of coupling, two branches 
of the dispersion curve come very close to each other. 
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Fic. 2. Phase velocity (h/a = 1/30, m = 2). 
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Fic. 3. Phase velocity (h/a = 1/30, m = 5 
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One such location is between the third and fourth 
branches (counted in the ascending order) of the solid 
curve. More detailed calculations were carried out for 
the case h/a = 1/30 and m = 2 around this location 
(66 < » < 76), and the results were plotted in Fig. 6 
with a much larger scale than that used in Fig. 2. In 
Fig. 6 we see that the line segments A;B; and C2D» are 
not collinear but approach a common straight line be- 
tween them as an asymptote, and the same is true with 
A.B, and C,D,;. A similar situation exists between the 
first and second branches of the dotted dispersion curve. 
Detailed calculations were made for the case h/a = 
1/30 and m = 2 for 56 < \ < 68, with the result plotted 
in Fig. 7. 

In contrast to Figs. 6 and 7, Fig. 8 shows quite a dif- 
ferent situation developed between the second and third 
branches of the dotted dispersion curve of Fig. 2. Re- 
sults given by the cubic frequency equation were rather 
irregular for values of \ between 103.6 and 104.2 in this 
case of h/a = 1/30 and m = 2. Indeed, within this 
range, numerical calculations yielded only one real root 
of 2 and thus also one real value of T for each value of 
\. The other two roots of 2 were found to be complex 
conjugates of each other, although their pure imaginary 
parts are much smaller than the pure real parts.* Con- 
sequently, there is only one, the lowest, branch of the 
dispersion curve within this range of \, and the other 
two branches as shown in Fig. 8 have a definite tend- 
ency of meeting each other both near points B and near 
points C, thus leaving the space between B and C un- 
connected by either of these two branches. Asymp- 
totes still exist, but they cover up the two branches 
from the top and bottom. This unexpected behavior 
of the dispersion curve according to the classical theory 
is probably due to the small k-terms neglected in Eq. 
(4.10) or due to rounding-off numbers in the machine 
calculations, which may have caused inaccuracy in the 
coefficients of the various terms and thus inaccuracy 
in the roots of the frequency equation. 


Frequency 
Much of the above discussion on VT also applies 


to VQ. For instance, in Fig. 9 is shown the result of 

* Complex roots of 2 have also been found in the other cases. 
For instance, the roots are 10862.330 + 46.124831 7 in the case 
of h/a = 1/30 and m = 5 when A = 104, and 1244.6020 + 


9. 


7.7233410 7 in the case of h/a = 1/10 and m = 2 when A = 35. 
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V 2 plotted versus \ for the case of h/a = 1/30andm = 
2. Since VQ is equal to V T multiplied by A, the 
Figure is virtually a distorted version of Fig. 2 obtained 
by tilting the curves in the latter higher and higher 
toward the right as A is increased, and many character- 
istics of the VI curves are seen te have been kept. An 
additional characteristic shared by all branches of the 


/ ° ° . . 
VQ curves is that they increase monotonically with \ 


Frequencies of infinite shells (and simply supported 
finite shells) were also calculated by Arnold and War- 
burton!! and by Baron and Bleich'* on the basis of 
equations excluding SR. The results were found to be 
in good agreement with those obtained here, but only 
a very limited range of values of \ was covered by these 
authors. Thus, in references 11 and 12, the values of 
do not exceed 4 and 7, respectively, and much of the 
interesting part of the frequency curve in the classical 
theory was not revealed. On the other hand, the con- 
clusion may be drawn that, within these limited ranges 
of X, the results given in these references are satisfactory 
even from the standpoint of the theory including SR. 


Group Velocity 

A typical graph of the group velocity is given for the 
case of h/a = 1/30 and m = 2 in Fig. 10, in which the 
various branches of the solid and dotted curves are 
identified by numbers in small solid and dotted circles, 
respectively. Previous comments on the range of 
agreement between the solid and dotted curves of the 
phase velocity, \ < 15 in the present case, may clearly 
be extended to apply to the group velocity. The limit- 
ing values of C/co for very small and very large values 
of \ also agree with those given earlier for the general 
case, but fluctuations are seen to exist for intermediate 
values of }. While the fluctuation is limited between 
0 and 1 according to the solid curve, it is much more 
severe in the case of the dotted curve. In fact, the 
group velocity not only increases with ) indefinitely 
according to the third dotted branch, but it also tends 
to become very large in both the positive and negative 
directions in the neighborhood of \ = 104 according to 
the second as well as the third dotted branch. This 
and another situation connected with the dotted group- 
velocity curve near \ = 61 are easily seen to be con- 
sistent with the findings shown in Figs. 7 and 8 about 
the phase velocity, if we keep in mind the relation be- 
tween the phase-velocity curve and the frequency curve 
discussed before and the relation between the fre- 
quency curve and the group-velocity curve indicated 
by Eq. (5.8). It is of course difficult to see how infi- 
nitely short waves may be propagated with an infinitely 
large group velocity, as was mentioned before. 


The group velocity for cylindrical shells does not 
seem to have been discussed before in the literature. 


Modes 


Graphs of the modes for the case h/a = 1/30 and 


m = 2 are shown in Figs. 11 and 12. The ordinate 
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U/W or V/W was plotted according to a method due 
to Donnell,!* which has been found to be best suited 
for our purpose. Thus, corresponding to the center 
region where L’/W and V’/W are nearly zero, the motion 
is dominantly in the radial direction. On the other 
hand, for regions near the top and bottom of the Figures 
where W/U or W/V is very small, the corresponding 
motion is more longitudinal or circumferential. By 
the method used here, the ordinate also shows clearly 
the sign (positive or negative) of U/W and V/W. 

Only three of the five branches of the solid curve are 
shown in Figs. 11 and 12, their order being again indi- 
cated by numbers in small solid circles. Similarly, 
the corresponding branches of the dotted curves are 
assigned numbers in small dotted circles. The second 
and third dotted branches are seen to be discontinuous 
curves in both Figures, in which each discontinuity is 
indicated by a pair of identically drawn small dashes 
or dots. It is interesting to note that, seemingly, the 
discontinuities may be removed by ignoring the sign 
of U/W and V/W and plotting only their magnitudes. 

In the case of the modes, the first branch of the solid 
curve and that of the dotted curve appear to agree 
with each other very well up to \ = 50, a value much 
higher than the corresponding one, A = 15, for the 
other quantities considered before. On the other 
hand, the agreement between the respective second and 
third branches is good only when X is smaller than about 
tor 5. For values of \ below 7 the modes have also 
been calculated by Baron and Bleich'* using membrane 
equations, and their results for the case of h/a = 1/30 
and m = 2 are found to be in close agreement with 
those obtained here for the same range of \ by more 
complete equations. This includes the limiting case 


of \ = 0, for which the following general expressions 
may be deduced from membrane equations: 

v 

— = 0, 2, 0 

WW 

V 1 2— (1+ pu)m? 

cn , —m 
if! m 2(1 + uw)m 


For the other limiting case of \ = ©, all the solid 
branches in Figs. 11 and 12 are seen to approach the 
centerline where U/W = V/W = 0. This is also true 
for the third dotted branch but not for the other two 
dotted branches. 


(IV) VIBRATIONS OF FINITE SHELLS 


(8) FINITE SHELLS WITH SIMPLY SUPPORTED EDGES 


For a finite shell of length /, the form of the displace- 
ments assumed in Eqs. (5.1) for the infinite shell may 
still be used provided that the parameter \ is now taken 
to be 


\ = nra/l (8.1) 


n being then the number of longitudinal half-waves 
contained in the length of the shell. 
the theory including SR the boundary conditions satis- 


According to 


SPACE 
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fied by the displacements in Eqs. (5.1) in this case of a 
finite shell are 


The first three of these conditions are obvious, and the 
last two may be verified by substituting Eq. (5.1) in 
the expressions of J/, and N, in Eqs. (2.4). On the 
other hand, if the classical theory not including SR is 
adopted, the conditions satisfied by the expressions of 
u, v, and w from Eqs. (5.1) are 


wnu¢g= MM, = N, = 0 
in which, according to the classical theory, 


M, = —D[(0°w/Odx?) + u(0?w/ds?) + (u/a) (dv/Os) 
N, [Eh/(1 — w?)] [(Ou/Ox) + w(Ov/0s) — u(w/a)] 


The conditions in both Eqs. (8.2) and (8.3) simulate 
those at a simple support. In both cases the non- 
vanishing longitudinal displacement u at the edge is 
permitted. The change of slope y in the longitudinal 
direction is further permitted to take a nonzero value 
at the edge according to the theory including SR; it 
naturally does not enter the problem when the classical 
theory is adopted. Physically these conditions imply 
that the shell at the simply supported edge is restricted 
to remain circular but is free to move in the longitudinal 
direction, and no force or moment is acting at the edge 
in this direction. 

Results for the frequency and modes obtained before 
for infinite shells may therefore be applied immediately 
also to finite shells with simply supported edges, pro- 
vided only that \ is now interpreted according to 
Eq. (8.1). 


(9) FINITE SHELLS WITH OTHER TYPES OF EDGES 


Vibration problems of finite shells with edge condi- 
tions other than that of the simply supported shell may 
still be solved, at least approximately, but conveniently, 
by making use of Galerkin’s method together with 
Donnell-type equations. For simply supported shells 
the use of Galerkin’s method would yield the same 
results as the above method of direct substitution. 
The method will therefore be illustrated by considering 
first the case of a finite shell with clamped edges. 

The natural frequency will now be computed on the 
basis of Eq. (4.10), which does not include SR. How- 
ever, instead of the original equation, a modified form 
of the equation obtained by application of the operator 
V~‘ will be used. Similar modification in the static 
case was adopted by Batdorf,'* who found that only 
through the use of the modified equation would con- 
vergent results be obtained for the buckling problem of 
cylindrical shells with clamped edges. For the same 
reason the modification is also adopted here.* 

For the vibration problem we substitute into the 
modified form of Eq. (4.10) the following expression of 
Ww: 

* It should be noted that such a modification in the equation 
used does not affect the result for the simply supported case 
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w= Wsin wt 


where IV is afunction of x and s. The result is 
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The deflection form used by Batdorf!* for the buckling problem of a clamped shell is also used here: 


W = cos (ms/a) 3 an, cos [(n — 1)ax 
n=1 


in which, as before, m is the number of circumferential waves, a the radius, and / the length of the shell. 


1] — cos [(m + 1)ax 


t 9 
l { (9.2) 


Each term 


of the series in Eq. (9.2) satisfies the following boundary conditions for w at the clamped edges: 


ow/0x = 0 (x = 0,x =)) 


a7) = 
vy = 


However, the governing equation (9.1) is not satisfied; when IV is substituted into the equation, its left-hand member 


does not vanish, and the equation may be put in the form F(x, s) = 0. 


the coefficients a, be determined by the equations 


omeiae wm ms (n — 1)xrx (n + 1)arx 
F(x, s) cos cos — cos dx ds 
0 /7 0 a l l 


The method of Galerkin then requires that 


[3 . See 


Carrying out the integration, we obtain the following infinite set of simultaneous equations: 


a,(2Mo + Me) — a3Meo = 0 
ao( My, + Mz) — aaM3 = 0 
Qn 0M, —_ 


where 

; , I ; ai I 
M, 034 ~—Fig? + -O~ sik” + k bias bho Pete, a4 
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5 (1 — w) + 5k] — wn ~~ =e A*¥ —k 

7-3 oe, Lo 
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These equations are homogeneous in a@,. To have nonzero values of @,, 


vanish; thus, 


a, (iM, 1 + Mas) + OntoM nasi = ( (n=. 


= U (nm = ay 
(9.3) 
54. 
— M)Y + 
—= - 
> : ee Ri — M)Y, xX 
}— 3m°r* — (4 — w)m'r* — m' (9.4) 


the determinant of their coefficients must 


ay, Qe as -_ ay a4 a6 
n= 1 2My + Mz —M, 0 = 0 0 0 
n=3 —M, M.+ M, —M, 0 0 0 
e=5 0 —M, M,+ Me 0 0 0 
: : ’ d = Q (9.5) 
sa Z 0 0 0 Ww4+ iM — IM. 0 
n= 4 0 0 0 —M; M,; + Ms; —M, 
n= 6 0 0 0 0 —M; Ms, + M, 
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The two subdeterminants may each be set equal to 
zero, leading to results corresponding to symmetrical 
and antisymmetrical modes, respectively. First ap- 
proximations may be obtained by keeping only one term 
in each of the subdeterminants. For instance, if only 


the a, term is kept in the first subdeterminant, we have 


2My + M2 = 0 (9.6) 


from which approximate values of the fundamental fre- 
quencies may be computed. The approximations may 
always be improved in accuracy by keeping more 
and more terms in the determinant. 

The form of u and v corresponding to that of WV chosen 
in Eq. (9.2) may be found by means of Eqs. (4.8) and 
(4.9). On the basis of the results thus obtained, the 
modes of vibration may be determined once the fre- 
quency is known. The boundary conditions on « and 
v implied by the method of solution used in this case 
have been discussed by Batdorf! for the corresponding 
buckling problem. 

The case of a flexibly supported shell may be treated 
similarly. The flexibly supported edge of a cylindrical 
shell refers to one at which the deflection w of the shell 
is zero and the slope in the longitudinal direction is 
directly proportional to the local bending moment .1/, 
in the same direction. When the simplifications adopted 
in reference 1 are introduced, the term (u/a) (Ov/Os) 
may be neglected in the following expression for the 
bending moment: 


M, = —D[(O°w/Ox?) + u(O°w/Os?) + (u/a) (Ov/0s) ] 


At the edge where w is zero, 0°w/Os? also vanishes, and 
M, becomes proportional to 0°w/dOx*. The boundary 
conditions of a flexibly supported cylindrical shell may 
therefore be written 
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in which the difference in sign in the two lines is due tg 
the usual sign convention used for J/,, and c is a pro- 
portional constant. The constant ¢ is an indication of 
the spring action of the support of the shell. 1, 
determine the value of c we must resort to approxi. 
mate means. For instance, in drum-type structures, 
the spring action of the end plate may be assumed to be 
such that it is determined by considering the case of 
axially symmetrical bending of a circular plate by 
uniformly distributed edge moment. Such being the 
case, it is easily shown that 


c = Da D,(1 + Ky) 


in which D, and yu, are, respectively, the flexural rigidity 
and Poisson’s ratio of the end plate. The value of ¢ 
may be determined for flange-type and other types 
of shell structures in a similar manner. It is clear that 
cylindrical shell structures can be more truly simulated 
by the flexibly supported case than by the two previous 
cases. Indeed, the flexibly supported case is much more 
general and reduces to the clamped and simply sup- 
ported cases when c takes the values zero and infinity, 
respectively. 

The following form of IV may be shown to satisfy 
the boundary conditions set forth in Eqs. (9.7): 
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Making use of Galerkin’s method on the basis of the 
modified form of Eq. (4.10) and following the same 
procedure as used in the clamped case, we arrive at 
the following infinite set of homogeneous equations for 
the coefficients a,,: 
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here /, is the same as given in Eq. (9.4). To have 
gonzero Values of a,, we must set the determinant of 
their coefficients equal to zero. This leads to the fre- 
wency equation, which, however, is much more com- 
oicated than Eq. (9.5) and will not be presented here. 
Ry keeping only the a, term in the series of W, the 
equation takes the following simple form: 


(4/3) (l/m*c)] + My [2(l/42c)? — 
+ M,[(1/4mrc)? — (1/3) (l/r’c)| = 0 


fl ai 


1/ xc) (9.8) 


om which approximate values of the fundamental 
irequencies may be calculated. Clearly, Eq. (9.8) 
reduces to Eqs. (9.6) and (6.1), respectively, for the 
damped and simply supported cases, on putting « 
equal to zero and infinity. 

It should be brought to the reader’s attention that 
in this treatment of the flexibly supported cylindrical 
shell the supporting system (for instance, the end 
plates in a drum structure) has been assumed to be 
massless, which will certainly impose restrictions to the 
:pplication of the result obtained. 

If desired, the more complete equation (4.6) may be 
used instead of Eq. (4.10) to determine the frequency 
of shells with clamped or flexibly supported edges. 
The procedure of solving the equation by the Galerkin 
method would, of course, remain the same. 

Numerical calculations for finite shells with clamped 
ind flexibly supported edges are expected to be carried 


out later. 
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Thermal Stresses in Plates 


M. J. Forray 


Principal Engineer, Structures Department, Republic Aviation Corp 
Farmingdale, N.Y 


June 16, 1958 


Sv” BULKHEADS are plate-like major components in the semi 
monocoque type of construction of air and space vehicles 
The temperature distribution is often unsymmetrical for vehicle: 


in various flight regimes so that 


temperature distribution exists on the skin and therefore on th 


bulkhead of the craft. 


It is the purpose of this note to determine the stresses in 


traction-free bulkhead subjected to a general two-dimensiona 


temperature distribution. A particular example of a 


bulkhead subjected to an approximate temperature distribution 
when a vehicle is flying through the atmosphere with a positiv 


angle of attack is, considered in detail. 


According to the two-dimensional, plane stress theory of elas 


ticity,! 


of (1) equilibrium, and (2) compatibility, at all points of the bod) 


plus appropriate boundary conditions 


strain, and displacement is found which satisfies all the required 
conditions, the solution is correct and unique! according to the 


uniqueness theorem of Kirchhoff 
In the plane-stress problem for simply connected domains onl 


one equation of compatibility! is considered—namely, 


(O7€r2/Oy?) + (O7€y,/Ox?) = O7%e7,/OxXOY 
If the stress-strain laws! 
€rr = (1/E) (err — voy) + aT (x, y | 
€yy = (1/E) (eyy — vorz) + aT (x, y 2 
éry = (1/G)ory, (O22 = O2r = Ory - 0) 


a rotationally unsymmetrica 


circular 


the components of stress must satisfy certain equations 


Then if a state of stress 
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ee , a, T represent the components of strain, stress, 


linear expansion, and temperature, respec 


and the equilibrium equa- 


coeflicient of 


re substituted into Eq. (1 


ns are used, there results the equation 

V%orr + Gyy) + EaV?T = 0 3) 
er V2 = (02/dx?) + (02/dy? 
This is the compatibility equation in terms of stresses and 


only If is introduced 


it the components of a stress are derivable from @ by 


emperature a stress function @ (x, y) 


o = 0°p/Oyv", Gy = 07/Ox*, ory = —O0*/OxXOY (4 
then equilibrium is satisfied automatically, and the compatibility 
equatic m on @ becomes 


Vig = —Ea V?T (o 


V4 = [(0?/dx?) + (0?/dy?)] [(0?/Ox?) + (0?/dy?) 


where 
The boundary conditions to be satisfied are that the tractions be 
zero at all points of the boundary curve of the simply connected 
jomain It is interesting and helpful to know that the bound- 
ury conditions of free traction can be expressed? in terms of the 


stress function @ which are 
o@ = 0¢/On = VU (6 


where 7 indicates the normal direction to the bounding curve 
A formal 
particular problem of a 


must be solved subject to conditions (6 


solution will 
circular plate subjected to a specific temperature distribution will 


Eq. (5 
now be given, and a 
follow 

The complete solution of Vig = 
To this we must add any particular solution of the 


0 is well known.! Denote it 


by oe = ¢@ 
nonhomogeneous Eq. (5). A particular solution of 


Vo => Eal é 
is necessarily a particular solution of Eq. (5 Assume that 
~EaT is expressible in a series of the form Por) + 35 [Palr 

n=1 
cos 00 Q,(r) sin mg] where P,(r) and Q,(r) are easily deter 


mined in any particular case and are assumed known. It is 














then readily verified that the particular solution of Eq. (7) is of 
the form 
oP = Sa) > g,(r) cos nO + h,(r) sin nO 8 
n l 
where 
gare f r-t Lf P,(r)ri-m dr dr, (n = 0,1, 2, 
. ~ 
Jy = » 2 ) 1 = 9 
} JS y JS Url’) ? drdr, (nm ye A 
Therefore @ = ¢ + Pp where 
Oo = ( Cuff? -- (Ciel > Cisl cos 6 + 
C1o'r + €13’r?)sin 8 4 BS (cyr” + d,r™*?) cos nO + 
n=2 
(c,’r® + d,,’r 2)sin né (9 


Terms are omitted from the ful] solution! in order to satisfy the 
condition that the stresses remain finite as 7 > 0) and @p is found 


from Eq. (8) 


CIRCULAR BULKHEAD 


Consider a circular bulkhead of radius a subjected to a temper- 
iture distribution 
T= 


T(r/a)? (1 —cos@) + 7; (10) 


where 7), 7; are constants and the polar axis (@ = 0) 
Thus 


is assumed 
to be vertically upward. 
Po = —EalT; + T(r a)?], 

P,; = EaT((r?/a?), Pr =0, n 


>1, @2.=0, x= 


Substitution of the expressions for P,,, Q, into Eq. (8) yields the 


particular solutions 


go = —(Ea/4) [T; r? + (Tor*/4 a®)) | 
—_ 11 
o, = Eal r™ 15a? { 


The combination of Eqs. (11) with Eq. (9) and satisfaction of the 


conditions on r = a results in 
QoQ = Eal ; (a?/16 e r*/8 r*/16a* 
ar/30 , 10a 4 r4/15a2)|cos 6 12 


The stress components corresponding to this stress function are 


Orr = EaTy\(1/4) [1 — (r?/a? + 
(r/5a r/a 1] cos 6} 
org = EaT((r/da) |(r/a 1} sin 6 13 
oo = EaTy j(1/4 l (3r7/a t 
(r/5a) {(4r/a 3] cos 6} 
and this is easily verified to be the correct solution 
REFERENCES 
! Timoshenko, S., and Goodier, J. N., Theory -f Elasticity, pp. 21-27 
55-58, 236-238, 421-425; McGraw-Hill Book Company, Inc., N.Y., 1951 
2 Boley, B., and Weiner, J., Thermal Stresses for Aircraft Structu WADC 
Technical Note 56-102, August, 1955 


Viscous Dissipation in Low Prandtl Number 
Boundary-Layer Flow 


E. M. Sparrow and J. L. Gregg 
Lewis Flight Propulsion Laboratory, NACA, C 
June 16, 1958 


leveland, Ohio 


IT A RECENT ISSUE of the JOURNAL, a set of approximate (but 


mathematically elegant) solutions of the energy equation for 
low Prandtl! Number boundary-layer flows was given by Morgan, 
Pipkin, and Warner The basis of their method lies in the fact 
that the thermal boundary layer is much thicker than the viscous 
boundary layer when the Prandtl Number is small. So, in the 
outer part of the thermal boundary layer, the velocity was taken 
from the potential-flow solution, and only in the inner part was 
an approximate correction made for the nonuniformity of the 
velocity distribution 

One of the problems studied by Morgan et al 
mination of the adiabatic wall temperature for dissipative flow 


over a flat plate the 


was the deter 


They found that excess of the adiabatic 


wall temperature, / over the free-stream static temperature, 
t_, is given as follows: 
(t teo)/(Ue?/2c,) = Pr'/? (0.9242 + 0.194 Pr'/2 1 


It is the purpose of this brief note to present supplementary infor 
mation on the adiabatic-wall problem in the form of exact (nu 
merical) solutions of the boundary-layer energy equation and of 
approximate solutions based on the Karman-Pohlhausen method 


The configuration is shown on Fig. 1 
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Fic. 1 
TABLE | 
Adiabatic Wall Temperatures 
bag = £1 T)* sae) 
Morgan, 
Pr Exact [Eq. (1) | 
0.03 0.1653 0.1659 
0.01 0.09434 0.09436 
0.006 0.07279 0.07272 
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KARMAN - POHLHAUSEN METHOD 
THIRD DEGREE POLYNOMIAL 


.2) FIFTH DEGREE POLYNOMIAL = 
-— ~ aa 
= ail 
OS 
aerT 
wait. 

Sow"! Pr Pl 

2 .08 r ee af Z 
Ug / 2c, —<EXACT SOLUTIONS 


a wet (NUMERICAL) 


MORGAN'S VELOCITY 


04 i APPROXIMATION, EQ (1) 
NACA 
02 
002 004 .006 Ol 02 04 
PRANDTL NUMBER, cpp /k 


Fic. 2. Adiabatic wall temperatures for flow over a flat plate. 


EXACT SOLUTIONS 


A similarity-type formulation of the adiabatic-wall problem 
for dissipative boundary-layer flow over a flat plate was first 
given by Pohlhausen in 1921. His analysis is described in detail 
by both Jakob? and Schlichting’ and need not be repeated here 
Numerical solutions of Pohlhausen’s energy equation have been 
carried out for Prandtl Numbers of 0.03, 0.01, and 0.006 on an 
IBM 650 computer. Results for the adiabatic-wall temperature 
are listed in Table 1 and are plotted in Fig. 2. (The abscissa 
range corresponds to the Prandtl Number range of liquid metals. ) 
Temperature distributions across the boundary layer are plotted 
on Fig. 3, on which the velocity distribution is also shown. The 
relatively large thickness of the thermal boundary layer is strik- 
ingly displayed here, verifying the underlying philosophy of 
Morgan’s method. The temperature distributions display the 
characteristic inflection-point profile of the adiabatic-wall prob- 
lem. 

For the purpose of comparison, Eq. (1) has been evaluated at 
the Prandtl Nuimbers of Table 1 and the results listed there 
The agreement between the approximate and exact solutions is 
seen to be excellent, the deviations being no more than a fraction 
of a per cent. 

To provide perspective for this noteworthy success of Morgan’s 
approximate method, it is logical to inquire as to whether his 
formulation leads to equally correct results for the situation where 
there is heat transfer taking place. A summary of low Prandtl 
Number heat-transfer results (without viscous dissipation) is 
given in reference 4. It was found there that for the case of uni- 
form wall temperature the Nusselt Number given by Morgan’s 
method is 3.7 per cent below that of the exact solution for Pr = 
0.03, with smaller deviations occurring at lower Prandtl Numbers. 
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Fic. 3. Temperature and velocity distributions for adiabatic 
wall situation. 
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For the uniform wall heat flux case, the deviations were found t 
be larger (and of different sign). 

It is interesting to examine why Morgan's method is not quit 
as successful in heat-transfer situations as in the adiabatic-walj 
problem. The first explanation which suggests itself is that th, 
thermal boundary layer may not be as thick in the heat-transfe; 
situation as in the adiabatic-wall problem. Inspection of avail 
able exact solutions shows that this is indeed true for the uniforp 
wall heat flux case. But, for the uniform wall temperature case 
the thermal boundary-layer thickness is not very different fror 
that of the adiabatic-wall problem. So, it would appear that 
difference in thermal boundary-layer thickness may not be the 
complete explanation for the uneven success of Morgan’s method 
Study of his paper shows that the results for the adiabatic-wal] 
situation lean more heavily on exact velocity solutions than d 
those for the heat-transfer situation. In the latter, only the 
asymptotic value of the Blasius function is taken from exact 
boundary-layer theory; but in the former, the values of two addi 
tional velocity parameters are taken from the exact theory 


K ARMAN-POHLHAUSEN METHOD 


It has been shown? that the Karman-Poiilhausen method ap- 
plies with reasonably good success to the heat-transfer problem 
for low Prandtl Number boundary-layer flows (Nusselt Number 
predictions correct to about 5 per cent Now, we turn to the 
application of this method to the adiabatic-wall problem 

A polynomial of third degree was first selected to approximate 
the temperature distribution. This profile, introduced into the 
integrated form of the energy equation, led to adiabatic-wall 
temperature results as shown by the lower dashed curve of Fig. 2 
The curve corresponds to a cubic velocity profile, but a recalcula- 
tion using a fifth-degree velocity profile gave essentially the same 
results. At Pr = 0.03, the Karman-Pohlhausen prediction 
falls about 16 per cent above the exact solution. This deviation 
increases with decreasing Prandtl Number to 38 per cent at Pr = 
0.008. 

With the idea of improving the performance of the Karman- 
Pohlhausen method, the temperature distribution was written as 
a polynomial of degree five. The result of this calculation appears 
as the upper dashed curve of Fig. 2. It is interesting (but not 
completely without precedent) that poorer results are obtained 
as the approximating polynomial becomes more complex 

We are forced to regard the findings of Fig. 2 as indicating one 
of the rare failures of the Karman-Pohlhausen method 


QUANTITATIVE EFFECTS OF VISCOUS DISSIPATION 


It is interesting to look into the practical importance of viscous 
dissipation on low Prandtl Number boundary-layer flows.  In- 
spection of Eq. (1) or Fig. 2 shows that effects of viscous dissipa- 
tion depend on only two fluid properties: Pr and c,. Among 
the technically important liquid metals (as listed in the Liquid 
Metals Handbook), the combination leading to the maximum 
effect of viscous dissipation is (approximately): Pr = 0.03, ¢, = 
0.03 B.t.u./lb.-°F. Utilizing these values in conjunction with 
Table 1, we are able to write that for any liquid metal 


lew ~— ta & 107* UVn?* (2 
~ isin ft./sec. and¢isin °F. Asa simple example, sup- 
pose that U. = 50 ft./sec., then tar — fo < 0.25°F. It thus 
appears that for low Prandtl Number fluids, the adiabatic-wall 


where L’ 


temperature problem is primarily of theoretical interest. 
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Summary of Low-Prandtl-Number 
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Flow Properties Behind Strong Shock Waves 
in Nitrogen? 


- F. Waldron E 
stitute of Aerophysics, University of Toronto, Toronto, Canada 
yne 16, 1958 


SYMBOLS 


sound speed 
Vu flow 
tatic pressure in mm. Hy 


Mach Number 


particle velocity 
static temperature 
1 hock velocity 
ibseript (1) refers to the state ahead of the shock wave; and (2) to that 


hind the way 


Dimensionless Ratios 


u 


HOCK WAVES were produced in a 2- by 2-in. shock tube by 
burning a mixture of oxygen, hydrogen, and helium at con 
stant volume The particle velocity, static pressure, and sound 


+ The work was done under the direction of Dr. G. N. Patterson and super 
Thanks are due to Dr. J. G. Hall, Dr. J. H. de 
Taylor, D. W. Boyer, and to the Computation Center, Uni 
The project was supported by a 


sed by Dr. I. I. Glass 
eeuw, B. W 
rsity of Toronto for their assistance 
rant from the Defence Research Board of Canada 

ow Head, Aerodynamics Section, CARDE, Valcartier, Canada 
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Fig. 1. Particle velocity (Ux) as a function of incident shock 
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Fic. 2. Variation of the pressure ratio (2) across the shock 
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speed were measured simultaneously with the shock-wave ve- 
locity for shock-wave Mach Numbers from 7 to 13 for an initial 
pressure of 25 mm. Hg of nitrogen 

The particle velocity was obtained from a schlieren wave speed 
photograph of the contact front produced by the interaction of 
the incident shock wave with a perforated plate The static 
pressure was measured with an S.L.M. quartz crystal gauge 
From the inclination of the schlieren traces of a spherical wave 
formed by the discharge of a weak spark, the sound speed was 
obtained. The spark discharge occurred in the region between 
the incident shock wave and the contact front. All measure 
ments were obtained at a distance of 15 ft. from the diaphragm 
nitrogen from the combustible 


which initially separated the 


gases 

The experimental quantities are compared in Figs. 1, 2, and 3 
with the theoretical flow quantities calculated in reference | 
where dissociation and vibration are considered to cause the 
deviation from the ideal (y = 1.4) properties. The measured 
values of particle velocity (Fig. 1) and pressure ratio (Fig. 2 
indicate good agreement with theory. From the measured values 
of shock wave and particle velocities, additional values of the 
pressure ratio were calculated by means of the equations of con 
tinuity and momentum, and these are shown on Fig. 2 
The flow Mach Number (Fig. 3 


measured sound speed and particle velocity 


was determined from the 
Additional points 
were obtained from the experimental values of particle and 
shock-wave velocities together with the theoretical values of the 
effective isentropic index rhe calculated values of the temper 
were obtained from the measured quantities 
Although a fair 


iture ratio (Fig. 4 
ind the theoretical degree of dissociation 
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amount of scatter exists in Figs. 3 and 4 the calculated values 
bracket the theoretical curves. 

A discussion of the results and further details are given in refer- 


ence 1. 
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Minimum Runway Length for an STOL 
Aircraft 


E. M. Shoemaker 
Boeing Scientific Research Laboratories, Boeing Airplane Co., 
Seattle, Wash. 


June 18, 1958 


SYMBOLS 


A = wing area 

CL coefficient of lift 

Cp = coefficient of drag 

D drag 

I (W/z) V(d0/dt) inertial force (Fig. 1) 
L = litt 

R (W/A) (2/pgC_z) radius of circular flight path 
7 thrust 

V = speed 

W weight 

h = altitude at point on flight path 

l = range after take-off 

7] = flight angle (Fig. 1) 

04 = maximum possible @ 

a thrust vector angle (Fig. 1) 

ri) air density 


W' WISH to examine the optimum climb capability for the 
class of STOL aircraft (T/W < 1) equipped with a vari 
able (rotating) thrust vector device. The assumption is made 
that the optimum aircraft climbs at constant speed and therefore 
constant lift and drag. Constant speed is possible because of the 


variable thrust vector device. This assumption is made not 





OPTIMUM _ S.T.O.L. 


| CONDITIONS: 
| T/W = | 
\ D=0 











Fic. 1. Circular flight path for optimum STOL aircraft 


SPACE 


SCIENCES NOVEMBER, 1958 
only for simplicity but because it is a good approximation to + 
climb conditions of practical STOL aircraft 


Referring to Fig. 1, the equations of motion after take-off ; 
the most general case of constant speed are 


D 4 
L+T7sina = Weos@ + (W/g)V(d6/dt) 


W sin @ = T cosa 


Eliminating a, the range and altitude measured from the ta] 
off point are then 


6; 
[= f V cos 6(dt/da)d@a = 
) 
"2 ¢) f a1 
0 


eA; 
h= J V sin 0(dt/d@)d6 = 
0 


cos 6 d0/F(6, L, D, 7 


(V?/g) f ~ [sin 6 d0/F(0, L, D, 7 
0 
where 
POL, D, T) (L/W cos 6 
V(T/W)? D/W + sin @)? 


Note that from Eq. (5) the constant speed condition can be mai 
tained only until 6 increases to the maximum value 


a D)/W { 


Oy = sin 


after which some change must be made in the flight path and/or 


speed. By comparing Eq. (1) and Eq. (6) it can be seen that 


6 = 04 corresponds toa = (0. Hence, upon take-off, a is a maxi 
mum and decreases to zero at 6 = 9y,. 

It is convenient to express L and D in Eq. (5) in terms of 
L = (1/2)pV?AC.,, D = 


(1/2)pV?2A Cp, where C, and Cp may be considered constants over 


V through the usual relations 


the low speed range to be considered. To obtain the minimum 
runway distance (ground run plus range /) for a given altitude 
h, realizable for the entire class of STOL aircraft restricted t 
constant speed climb, we consider the optimum aircraft wheré 
thrust is maximized (7'/W = 1) and drag is minimized (Cp = () 
Eq. (5) becomes F = pV?AC,/2W and Eqs. (3) and (4 


l = (2/pgC_) (W/A) sin 6 i 


cos 6 N 


h =(2 peC,) (W A) {1 


so that the flight path is independent of V and take-off is possible 
at infinitesimal speed. If now V — 0, the ground run will ap 
proach zero and the total runway length is equal to the range 
It will be pointed out that these results are not obtained if 7°/W = 
land V — 0 without first setting Cp equal to zero 

Note that 05, = 90°, and that the flight path is the arc of a 
At the top of the 
are the aircraft would be capable of continued vertical flight as 
Also, Z = J, and, from Eq. (1) with D = 0 


+ @ = 90° so that the thrust is always vertical 


quarter circle of radius R = (2/pgC,) (W/A) 


indicated in Fig. 1. 
and T = W,a 


The range to any given altitude is from Eqs. (7) and (8), 


L= ~V/2hR —h?, h<R (9 


which represents the shortest runway length to clear a given 
obstacle height for the class of STOL aircraft 
ately from Eq. (9) that for any STOL aircraft the runway dis- 


tance must always exceed the obstacle height. 


It follows immedi- 


In Table 1, radii R are tabulated for various combinations of 
lift coefficients and wing loadings in the region of practical in 


TABLE | 
Radii R of Flight Paths and Runway Length / Required to Clear 
a 50-ft. Obstacle for Various Values of C, and W/A 


R(ft 1(ft — 
CL/ Ww A 100 50 25 100 50 25 
1.5 1750 873 437 415 291 203 
2 1310 655 328 358 251 174 
! 655 328 


164 251 174 117 
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t:; optimum runway lengths required to clear a 50-ft. ob 


rest, ! 


tacle, frequently regarded as a standard, are also presented. It 


in be concluded that runway distances of 200-300 ft. (obstacle 
ight of 50 ft.) which have been frequently mentioned as within 
e reac f practical STOL are not attainable except at great 
acrifice in the form of large power plants and large wings or 


ther high lift devices 

ind (4) may also be integrated under the assumptions 
land ’ ~0O. The results, which will not be pre 
ire complicated expressions involving both C;, and 


Egs. (3 
that 7/ Vi 
sented here, < 

and a class of flight paths for fixed C, and W/A and variable 

» of which the path corresponding to Cp = O is the optimum 
ircular path presented here. The results for nonzero Cp, while 
jot representing the optimum condition since the values in 
Table 1 are increased, are more representative of practical con 
jitions. The integration of Eqs. (3) and (4) for 7/W < 1, corre 
sponding to a practical STOL aircraft, has not been accomplished 
nclosed form 

Studies based on the results of the optimum flight path have 
proved useful in distinguishing between the competitive realms 


fSTOL and VTOL 


An Analytical Procedure for Orthogonalization 
of Experimentally Measured Modes 


Sidney |. Gravitz 

Senior Dynamics Engineer, Dynamic Science Section, 
American Aviation, Inc., Columbus, Ohio 

June 20, 1958 


North 


SYMBOLS 


Cij structural-flexibility influence coefficient expressing the de 
flection at station 7 due to a unit load applied at station 7 
A generalized stiffness influence coefficient expressing elastic 
coupling between modes 7 and 7 
generalized mass matrix element expressing inertia coupling 
between modes i and j 
physical mass assigned to station 7 
total number of mass control points on the physical structure 
number of measured modes 
6 measured mode shape 
orthogonal mode shape 
w modal frequency 
matrix 
diagonal matrix 
matrix inverse 
matrix transpose 
{ column matrix 
x matrix of order a rows by 6 column 


DYNAMICIST makes use of many analytical procedures 


her 
which require orthogonal mode shapes and frequencies as in 











put. One primary advantage of orthogonal-mode usage is that 
the accuracy of the theoretically determined elastic character- 
istics can be experimentally checked during ground vibration 
tests, and, in fact, the experimentally determined modes should 
be directly usable in subsequent theoretical analyses 

Unfortunately, it is well-known that experimentally determined 
ground vibration test modes are invaribably not orthogonal. The 
following questions immediately arise: 

a) How nearly orthogonal are the measured mode shapes? 

b) How may the measured modes be modified or converted 
into corresponding orthogonal ones? 

The procedure described herein utilizes measured vibration 
data plus corresponding mass data (usually theoretical) to com 
pute a matrix of structural flexibility influence coefficients. The 
closeness of this matrix to being symmetrical is a measure of the 
closeness of the input modes to being orthogonal. The pseudo- 
experimental influence coefficient matrix is then made symmet- 
nical by taking the average of the off-diagonal elements. This 
Averaging process is equivalent to that performed when reducing 
data from an experimental static influence coefficient test. The 
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symmetrical influence coefficient matrix plus the input mass data 
are then used in a conventional vibration analysis to obtain 
orthogonal modes. Thus the orthogonalization procedure may 
be considered to perform two functions 

(a) Averaging out of imperfect experimental techniques, in- 
strumentation errors, and structural nonlinearities 

(b) Averaging out of nonlinearities introduced by large meas- 
ured amplitudes being used in theoretical procedures which as 


sume very small amplitudes 
DESCRIPTION OF PROCEDURE 
The following is assumed given: 
m x es 9x “ x 


Starting with the classic free-vibration equation 


io} = w2(C m ot 
or, summarizing the equations for a set of ‘‘g’’ modes 
olanxg = [Clax mM~|nxnl@jn xq) ~w*~]¢ xq 


Solving for the matrix product of the influence coefficients and 
the masses: 

c m = |¢| w? o 
influence coefficients can be 


At this point it is noted that the 


straightforwardly obtained as an ‘‘n x n’’ matrix, if ‘‘n’’ modes 


are measured at ‘‘n’’ control points. However, for most practical 
cases, it is desired to measure relatively few mode shapes at rela 
tively many control points (m >q This results in a rectangular 
mode shape matrix with an inverse mathematically undefined 


To circumvent this difficulty, it is here assumed only temporarily 


that m = q so that the matrix inverse is mathematically defined 
In order to obtain an expression for the mode shape matrix in 
verse, we make use of the generalized mass 
[Wn ol" m o 
Solving for the desired quantity: 
lp = [M]—! [o]° [ym 


Substituting back into the previous equation 


[C] [-m ] = I¢ 2 (mM o n 


(\o m o xe 


It is now possible to itemize the steps in the procedure: 


(1) Compute 


K measured > mM w* 
7x4 = 1 ; “ 
= lox m x Ol, w ee 


(2) Compute 


C |measured = 
nx) 


3) Obtain 


. - 9 . 1 ’ 7 
C jmeasured = ~ C |measured C |measured 


symmetric 
uxn 
Averaging the off-diagonal terms is but one of several means 
a symmetrical influence coefficient matrix. This 
orrespondence 


of obtaining 
particular approach was chosen because of its close c 
with the experimental static influence coefficient determination 
procedure 

4) Substitute into the free-vibration equation and utilize such 
conventional procedures as matrix iteration or frequency choice 


to obtain orthogonal modes: 


io} = w?[C] [mn] 14} 
or ({C! [~m — (1/w?) [I]) (o} = 40} 
DISCUSSION 
It is noted that the procedure has certain reasonable limita 


tions: 
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(a) The number of meaningful output orthogonal modes can- 
not be expected to exceed the number of input modes (although 
in every case investigated to date, the number of good output 
The ‘‘exact”’ 


modes was equal to the number of input modes). 


influence coefficients incorporate contributions from the entire 
The pseudoexperi- 


infinity of modes for a continuous structure. 


mental coefficients are merely approximations to the 


ones, of accuracy sufficient to yield back data regarding only 


the modes that went into their formation. 


(b) If a set of forced responses bearing no relation to the true 


modes of the structural system is used as input, then the ‘‘orthog- 
onal’? output modes are physically meaningless. This is_be- 
cause the description of the physical system would be incomplete 
without specification of the imposed restraints. 


straints would in general be unknown from a conventional vibra- 


These re- 


tion test when data far from the frequency response peaks are 


used. 

(c) When considering a complex structure, such as a com- 
plete aircraft, and data concerning only one component (the 
wing for example) are desired, then sufficient control points must 


be included to adequately account for the generalized mass of 


the entire vibrating structure. In one practical case, it turned 
out to be adequate, when considering the wing, to account for the 
rest of the aircraft by means of only two lumped parameter con- 
trol points. 

The procedure described above has been coded as a production 
digital computer program. It has thus far been applied with 
success to several full scale aircraft components, as well as a 
complete aircraft. Typicaliy, the orthogonalized frequencies 
differ by no more than a few per cent from the input frequencies, 
and the orthogonalized mode shapes differ by less than 10 per 


cent at corresponding stations from the input shapes. 


A Note on the Interaction of a Normal Shock 
Wave With a Thermal Boundary Layer 


John P. Appleton and Hubert J. Davies 


Department of Aeronautical Engineering, University of 
Southampton, England 


July 9, 1958 


—— TWO-DIMENSIONAL FLOW considered is that of an ideal 
compressible fluid bounded by a rigid nonconducting wall of 
infinite extent parallel to the main-stream direction. A shock 
wave normal to the rigid boundary divides the flow into sub- 
sonic and supersonic regimes. In the supersonic region, a ther- 
mal boundary layer, in which the temperature distribution is in- 
variant in the direction of the mainstream, disturbs the normal 
shock wave and thus causes a modification to the pressure and 
velocity of the gas in the subsonic region. 
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The linearized theory contained in this note follows a treatmey 
developed by Lighthill,' and applied to the problem, which js ¢} 
subject of this note, by Griffith.2 In his paper Griffith shoy 
that the perturbations in pressure and velocity in the subsopj 
region are derivable from a potential function taken in the for; 
of a distribution of sources normal to the rigid boundary ap 
directly proportional to the excess temperature in the therm, 
layer. This note provides an alternative method which yie| 
analytic expressions for the pressure and velocity perturbation: 


THEORY 


The theory is presented in detail from that point at which tt 
method diverges from that of Griffith? in which report details 
notation and the preliminary work will be found. Reference j 
made to orthogonal x, y axes where the y-axis is taken to repr 
sent the position of the undisturbed shock and the x-axis ¢} 
rigid boundary. The two-dimensional steady flow equatio: 
expressing conservation of mass, momentum, and energy ar 
linearized in the usual manner. From these equations the per 
turbations in density and the component of velocity paralk 
to the x-axis are eliminated yielding two equations satisfied by th 
pressure perturbation p(xy) and the perturbation (xy) in tl 
component of velocity parallel to the y-axis. By defining a se 
of reduced coordinates these equations reduce to the Cauchy 


Riemann equations 


Op*/dx* = dv*/dy*, Op*/oy* = —(Ov*/dx*) 
Thus, a complex function W = p* + iv* may be defined which, if 
the inclination of the shock is ignored, is analytic in the domair 
x* > 0, y* > 0 of the z* (= x* + iy*) plane. Asv* = 0 onthe 
boundary y* = 0, the domain may be extended by reflection in th 
real axis so that W is analytic in the half plane x* > 0. The trans 
formation 

gc=(1—2*)(1 +: 2 


conformally maps the half plane x* > 0 onto the interior of the 
unit circle |¢| = 1. Onthe boundary ¢| = 1, the real part of 
W(¢) has the value p*(6} which is determined by the conditions in 
the supersonic region of the z* plane (see Fig. 1). In this region 
a general, sectionally invariant, temperature distribution of the 
thermal boundary layer is represented by 7(y) = 7\|1 + T’(y 
Ty) = 1. 
that, correct to the first order of the slope of the shock, p* = 
A T"(y), which is taken on x* = 0 


From the oblique-shock-wave relations one finds 


where A = —[2/(y + 1)] (Pi/P2) 21 — M.2)!/? 
Hence, as p* is periodic with respect to 6( = arg ¢), the condition 
on the unit circle |¢| = 1 may be written in the form 
p*(0) =A D> T,' cos n0 
n=O 
p*(@) = Rj A pe T,,’ fo" |, R = real part of 
n=0 

where ¢ is the value of fon |; = 1. It follows that 


We) = p* + iv* =A > 1,5 3 
n 0 
satisfies all the conditions of the problem. Thus, from Eqs. (2 
and (3) the pressure and normal velocity perturbations can be 
determined at any position in the subsonic region 
Following Griffith,? the profile of the disturbed shock is give 


approximately by 
By) = S &o dy = cz (0, y) dy 


where C = (U,; — U2)~!. By using Eq. (3), this may be written 


as 


B(@) = Df z. T,,’ sin n@- sec? (0/2)d6 4 
n 0 


where 
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Evaluation of Eq. (4) together with Eq. (2) will then yield an 
expression for the shock-wave profile B(y). 

The theoretical results thus obtained are compared graphically 
with theoretical and experimental results obtained by Griffith,? 
see Fig. 2. Without knowing in detail] the particular apparatus 
und under what conditions the experimental results were ob- 
tained, it is difficult to comment on the scatter of the experimental 
points around the theoretical shock-wave profile. 
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fa Shock Wave with a Thermal Bound 


Approximation of the Eigenvalues for Heat 
Transfer in Turbulent Tube Flow 


Robert D. Cess 

Heat Transfer Section, Mechanics Department, Westinghouse 
Research Laboratories, Pittsburgh, Pa 

June 20, 1958 


RF CONVECTIVE HEAT TRANSFER in laminar tube flow, Sellars, 
Tribus, and Klein! have obtained an asymptotic expression 
for the eigenvaiues by use of the WKB approximation. It is of 
interest to determine the possible application of this method to 
turbulent flow, and although more complete solutions are avail- 
ible,» * the analysis of Latzko‘ serves as a convenient example. 
Latzko has shown that under certain simplifying assumptions, 
which are approximately applicable for a Prandtl Number of 





FORUM 


unity, the energy equation yields 
(d/du) {((1 — u?) (dR/du + AVR = 0 (1 
with the boundary conditions 

R(O) = 0, R"(1 <M 
The variable u is defined in terms of the dimensionless radial dis- 
2\1/7 


tance as u = (1 — 7 ; the function R represents the radial 


temperature distribution; and A is the eigenvalue 
Following Sellars, Tribus, and Klein, a solution of the form 


is considered, where 


and, since A is assumed to be large, only the first two terms of the 


above series are retained. Upon substitution of Eqs. (2) and (3 


into Eq. (1), and equating powers of A 


l 
gy) = “i f V e/ £7) dt 


= —In{u’/4(1 11/4 


S 


such that R is given from Eq. (2) as 


1 y 
R= aexp [aif VY 7/1 eae] : 
l / = 
B exp | - vi f Vi/(l ay ae | u?/4*(] u7)'/4 (4 


Employing the change in variable y = 1 u, Eq. (1) becomes 
for small y 
(d/dy) |wWdR/dy)| + (A2?/7)R = 0 
and taking R = 1 for y = 0, the solution of this equation is 
R = Jy [2x V(y/7)] 
which may be written in its asymptotic form as 
R = (1/xd)!/2 (7/y)"/4 cos [20 -Wy/7 — (4/4 (5 
Correspondingly, for small y Eq. (4) becomes 
R= [4 exp (2\i Vy 7) +4 
Bexp (-24 V4 7)| (7y)'/4 (6 
and comparing Eqs. (5) and (6) the coefficients A and B are deter- 


mined such that Eq. (4) gives 


R = (7/md)!/? X 
e | a 
” V/A £7) dé (r/4 ] n/4(1 u7)! t 


= ‘ us 
\ 

cos | A 
‘ L 


For small «u Eq. (1) becomes 
(d?R/du?) + #uw?7R = 0 


for which the solution is 


R = Du}/? Jy /9[(24/9)u®/?) 4 Eu)/*J_1/9[(24/9)u®/*| (8 
This equation may be written in its asymptotic form as 
R = (9/ad)!/? (1/ut/*) {D cos [(24/9)u”* — (107/36)] 4 
E cos [(2/9)u9/* — (82/36)|} (9 
whereas for small « Eq. (7) becomes 
R = (7/md)'/2 (1/u7/4) cos [—(2d/9)u9/? 4 
KA — (7/4 (10 
i. 
where K = { fe/(l £7) dt 
. 
which may be expressed in terms of gamma functions as 
K = (WV 2/7) [1(9/14)/T(8/7)] (11 


Determining D and E from Eqs. (9) and (10), the solution of R 
valid for small u is found from Eq. (8) to be 
R = 2u!/2 sin (#/18) {sin [AK — (172/36)]Jiyo [(2A/9)u*’*] - 


sin [AK — (192/36)|J —~—1,/9|(24/9 u®/=|} (12 
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TABLE | 
n Eq. (13) Durfee 
0 1.38 2.9542 
l 2.68 12.346 
2 20.98 20). 858 
3 29.28 29 .252 
4 37 . 58 37.605 
5 15.87 15.94 
6 54.17 4.26 
7 62.5 62.6 
bd 70.8 70.9 


In order to satisfy the boundary condition R(O) = 0, Eq. (12) 
requires that 
AwK = [n + (19/36)] xr 


or rn = [7 V2 1(8/7)/T(9/14)] [n + (19/36)] (13) 


which is the asymptotic expression for large \. A comparison of 


this result with the eigenvalues computed by Durfee® is shown 
in Table 1. 

From this comparison it is evident that the method of Sellars, 
Tribus, and Klein may have definite application to determining 
the higher eigenvalues for heat transfer in turbulent tube flow. 
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On the Length of Hypersonic Nozzles 


James T. Kenney and Ying-Nien Yu 

President and Specialist, Respectively, Sandberg-Serrell Corp., 
Engineers, Pasadena, Calif. 

June 24,1958 


IT THE DESIGN of supersonic wind tunnels, it is customary to 
contour the nozzle wall until the flow is theoretically uniform 
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and parallel. The ordinates of the physical wall are obtained | 


adding to the perfect-fluid ordinates (7) the amount of the | 
boundary-layer displacement thickness(6* Thus, correct all 


ance is made for mass flow at every station. This method 
equally applicable to hypersonic nozzles and successful results 
reported. 

In hypersonic nozzles the length-to-exit-half-height (or radiy 
in axisymmetric case) ratios are much larger than in the super 
sonic case the Mach angle 


length combined with the generally high supply temperature r 


due to low Moreover, t gre 


sults in thicker 6* growth. Shortening the nozzle by omitting 
section near the exit will not only reduce this excessive boundar 
layer growth, but will also optimize the size of the usable 
region 

The usable test region lies inside of the actual boundary la 
The tabulated data in references | and 2 show that this Ja) 
Mach Number 


Hence the maximum usable test-section height is located at 


(6). 


range of 6 1? 


is approximately 1.56* for a 


point given by the following relation: 
d(r + 6* 6)/dz d(? 0.56*)/dz = 0 


Downstream of this point the test region actually begins to d 
crease in height 
A more conservative approach to test-section height opti 
ince 2(6 6;*) cot yz 
This is 


cteristic 


zation is to terminate the nozzle at a dist 


6,* cot w; upstream of the theoretical exit ipproximatel 


the point where the right-running char wave, whi 
meets in a common point with the exit characteristic wave a 
the inner surface of 6, originates from the perfect-fluid w 
The measurements on existing hypersonic wind tunnels indicat 
that 6,*/7,; can be as great as 1/4, therefore, the cutoff lengtl 
approximately (.1/\7/4 


Fig. 1 shows the cutoff length for a hypersonic nozzle. 
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The Practical Significance of Limit Analysis’ 


Philip G. Hodge, Jr 


Professor of Mechanics, Illinois Institute of Technology, Chicago, 


June 29, 1958 


ys THEORY of limit analysis is based upon the idealization 
known as a perfectly plastic material. Such a material hasa 
sharply defined yield stress under which the strains can increas¢ 
indefinitely. In particular, if the material is rigid perfectl; 
plastic, it is assumed that no straining is possible for stresses 
below the yield stress. Let 
tion of a set of forces determined to within a magnitude factor / 
As P is slowly increased, the unique value of P for which the rigid 


a given structure be under the ac- 


perfectly plastic structure can first undergo any deformation 1s 
variously known as the limit load, yield-point load, or collaps 


load 

If the structure is made of an elastic perfectly plastic material 
it can strain elastically for stresses below the yield stress and 
For such a material, the limit 


le- 


can 
never support any greater stress 
load is defined as that value of P for which indefinitely large ¢ 
formations could occur if geometry changes could be neglected 
It can be shown! that this definition leads to the same value as 
that for the rigid perfectly plastic structure 

Most real materials undergo a certain amount of strain hard 
ening and can, in fact, support stresses greater than the yield 


stress. Fig. 1 shows experimental stress-strain data for two struc 


* This investigation was sponsored by the Office of Naval Research 
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a turally important materials—mild steel? and 245-T aluminum 
30 Both materials have quite sharply defined stress levels at whi 
— the behavior deviates from the linear. In the case of aluminum, 
°8 ° ° 1 1 
“3 °) strain hardening begins immediately; for mild steel, there is a 
substantial increase in the strain before strain hardening. I) 
each case, the elastic behavior is given bv Hooke’s law, @ Ke, 
and the strain-hardening portion is approximated by a power 
function o = Be’ Table 1 shows the physical constants as 
ciated with the two materials 
rhe symmetric three-bar truss shown in Fig. 2 has been ana 
lyzed! for the two real materials and their corresponding elastic 
perfectly plastic materials; the results are shown in Figs. 3 and 4 
rhis truss has frequently been used in the literature,® iS pro 
viding a simple example which nevertheless contains many essen 
tial characteristics of more complex structures 
if The significance of the perfectly plastic limit load will be dis 
5 Experimental points. for mild steel cussed first in relation to the aluminum trus In the fully elas 
@ Expenmental points for 245-7 aluminum” tic range, an increase of 10 per cent, say, of the load is accompa 
- Fitted stress-strair arves “ ? : - 
. ‘ 7 nied by an increase of 10 per cent in the elongation \t d of 
Perfectiy plast stress-strain curves 
76,700 Ibs. the central bar of the truss becomes plastic Po ob 
i = ———— "/ tain a 10 per cent increase in load to 84,400 Ib e displacement 
02 20% 006 2.0 0.10 0,12 14 0.16 018 - : m - . : , 
must increase from 0.043 in. to 0.054 in., or 28 per cent. Thus th 
€ in r -— ° . ° . 1 
; ; ratio of displacement increase to load increase almost triple 
Fic. 1. Stress-strain curves for aluminum and mild steel , | , 
when the truss becomes partly plasti 
For a load of 107,900 Ibs., the remaining bars become plastic 
The displacement at this point is 0.082 in \ 10 per cent in 
PARLE | crease in the load now brings it to 118,700 Ibs. It follows from 
Physical Characteristics of Mild Steel and Aluminum Fig. 3 that the displacement increases to 0.164 in., or by 100 per 
Young’s modulus E 27.2 X 10° psi 12.51 X 108 psi cent Thus the displacement-rate load-rate ratio which less 
Yield stress o* 36.8 X 10% psi 45.0 X 10% psi than tripled in the partly plastic range now more than triples 
Strain at initial hard Bie a again and becomes ten times its elastic value 
ening € 0.0017 in. /in i ; k : 
; 4 Ps ot . se tO 1e steel truss, tl S Ss no mor iking n tl 
Hardening modulus B 125.4 X 10% psi 94.65 X& 10° psi For the steel truss, the results are eve © striking sane 
Hardening coefficient n 0.275 0.136 range of contained plastic deformation a 10 per cent increase in 
load requires a 25 per cent increase in displacement \t the in 
stant that P first reaches the limit load, 89,000 Ibs. the displace 
ment is 0.0324 in., a 10 per cent increase in load requires a dis 
H she H " placement of 0.268 in., or an increase of 827 per cent. Even if 
6 hd this latter displacement is compared with the displacement at the 
ae | . ; : 
7% | onset of hardening the increase is over 114 per cent 
\ 
In conclusion, the limit load which has a precise meaning for 
\ 3 | 2 the perfectly plastic material, has been shown to have qualitative 
i. , H significance for more realistic materials Below the limit load, 
eS r the deformation increase associated with a given load increment 
» | ° e ° ° 
* / is of the same order of magnitude as in elasticity, whereas above 
\ | A 
NY. the limit load it is greater by a factor of ten or more. It follows 
¥ that if the load on a structure is restricted to values less than the 
limit load, the deformations will remain of the elastic order 
\bove this load substantially larger deformations must be ex 
pected 
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Fic. 2. Symmetric truss ' Peer Wank eles. © Cte. There Perfectly Plastic Solids. 1 
Wiley and Sons, Inc., New York, 1951 
H —— 
& ——— an vor 
: r = 
: j a 3 
; - - ~------- : , 
4 
4 
/ — aluminum truss 
pertectiy plost russ " m wes 
[ *« - 
« O10 5 F O28 % ase a ro “Sie 
Clongeation nches 
Fic. 3. Load-elongation diagram for truss at Fig. 2, aluminum Fic. 4. Load-elongation diagram for truss at Fig. 2, mild steel 


= 12in., A; = A, = 1 in.’ H = 12in., A; = A 1 in.? 
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Or™= AN AEROELASTIC SOLUTION of a wing with a tip control 
surface is desired. For this case a noniterative superposi- 
tion method of analysis has been presented by Brown, Holtby, 
and Martin.'! The difficulties encountered in this type analysis 
are: (1) the solution of the flexibility matrix for the entire wing 
structure can get out of proportion to the capacities of digital 
computers and (2) the sectional properties between the control 
surface and the rest of the wing are sufficiently different to make 
the single cubic twist a poor approximation for the entire wing. 
A sectional aeroelastic solution which alleviates these difficulties 
is now described. 
Using the following angle-of-attack symbols: 


a; = initial angle-of-attack distribution 
apr = elastic angle-of-attack distribution 


ap = angle-of-attack distribution at which equilibrium is 


established 
then 
ar =ar+ar (1) 
Letting 
7 = nondimensional semispan wing station, (; = 1.0 at wing 
tip) 
ne = nondimensional semispan wing station of inboard end of 


tip control surface 


Two sections of the wing are established and the span station of 


each section is denoted by a subscript. Thus, 
m =7, When O0O<y< (2a) 
nm =n —n, When y»<n< 1.0 (2b) 
Sectional angles of attack are identified by superscripts 
a’ =a, when 0<7n<7” (3a) 
a” =a, when 7, 7 < 1.0 (3b) 
The curve-fit coefficients are defined as 
a = linear portion of elastic twist in region 1 
b = parabolic portion of elastic twist in region 1 
¢ = cubic portion of elastic twist in region 1 
d = constant portion of elastic twist in region 2 
eé = linear portion of elastic twist in region 2 
f = parabolic portion of elastic twist in region 2 
g = cubic portion of elastic twist in region 2 


Since the sectional properties of the wing and control surface 
are assumed known, the spanwise twist due to any loading can be 
computed and curve fitted in each section. For the initial wing- 
control surface loading the twist is 


a = dm + bon? + com? 


(4a) 


(4b) 


ay” = dy + eon: + fom? + gon? 
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By assuming equilibrium angle-of-attack loadings the elastj 


angle-of-attack distributions can be computed similar 


(4) 


then 


Similarly 


and 


and 


and 


and 


and 


The elastic angle-of-attack distribution expressed to a 


Letting 


ath = 7 
ag’ = am + bin? + cin 
ap” =d, + ein. + fine? + gin 
ar = 2 « 
are” = dym + bem? + com 
ag; = d; + e2n2 + fone? + Jone 
ar =7 
Qa} = 423m. FT b m + C31 
” 
QE = d Tt €3n2 TT S3m2* + 23 
ars” = 1.0 
an’ = ayn + byni? rT 4am” 
ap” =d,+ esn2 T ts 2? + guano 
” 
QPS = 4 
ars’ = asm + bem? + com? 

” , ‘ . 
aps = ds + esne + fone? + gsm 
ars” = ne? 
are = dem + bem? + cen’ 

” > «6 
ar = dg + ene + Jem" + Lene 
ar,” = 
agi’ = aim + bm? + cm 


fone? + gin* 


Si 


ap7” = d; + em 4 


in each section is 


ap’ = Am + Bn? + Cy 


ap” = D+ En + Fh? + Gm 


Substituting Eqs. (5), (7), and (9) in Eq. (19a), 


, 


QE a Car é 


= Aar;’ + Bape’ 


to Eqs 


(18a 
(18b 


cubic 


(19a 


(19b 


(20a 
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Substituting Eqs. (11), (13), (15), and (17) in Eq. (19b), 


apg” = Dar” + Ears” + Fare” + Gar:” (20b) 
Substituting Eqs. (20a) and (20b) in Eq. (1), 

ar’ = ay’ + Aap’ + Bar’ + Car;' (2la) 

ap” = ay” + Dar," + cars” + Fars” + Gar,” (21b) 


Now writing the initial angle-of-attack distribution for each of 


the final equilibrium angle-of-attack terms in Eqs. (2la) and 
21b), 
(1 — ay — de — a3 —a 

—b, (l — bs) —b; — bs 

—C) —C2 (1 — ¢3) —( 

—d —d» —d. (1 d 

—é —€» —¢ —¢ (1 

=o —Js Je fs 





The coefficients A, B, C, D, E, F, and G can be determined from 
Eq. (23) and the aeroelastic solution is determined by substi- 
tuting these coefficients in Eqs. (19a) and (19b). 

REFERENCE 
F., and Martin, H. C., A Superposition Method 
Behavior of Swept Wings, Journal of the Aero- 
8, pp. 531-542, August, 1951 


1 Brown, R. B., Holtby, K 
Calculating the Aeroelastic 


nautical Sciences, Vol. 18, No 
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A’ SCHLIEREN PICTURES taken in supersonic wind tunnels de- 
crease in quality when the density in the test section be 
comes low (lower than about 1 to 10 mm. Hg depending on the 
sensitivity of the schlieren setup) attempts have been made by 
several people to make the flow visible by self illumination. The 
most successful method so far proved to be the afterglow tech- 
nique developed by Kunkel and coworkers! which is especially 
useful for extremely low pressures (0.1 mm. Hg and lower). 
Another method which has received only little attention is the 
glow discharge method.? Recently, however, the author found 
some encouraging results with a modification of the latter method 
Though still of preliminary nature they indicate clearly that this 
particular glow-discharge method holds some advantages over 
The experimental setup is very simple 

A metal model which is held by an 


other similar methods 

and is sketched in Fig. 1 
electrically insulating sting carries a positive voltage of several 
thousand volts with respect to the grounded tunnel walls. These 
metal walls have been insulated around the model by a lucite in- 
sert which carries a thin steel frame at its front end to increase its 
mechanical strength and also to take up the discharge which 
otherwise could slowly deteriorate the highly polished surfaces 
of the wind-tunnel walls. The discharge is struck between the 
front end of the model and the steel frame and is blown down 
stream by the flow. Thereby, it spreads over a fairly wide region 
shows a photograph of such a discharge taken on Royal 
4.5 and 1 The current in this case was 
about 20ma. and the dissipated power about 40 The 
Mach Number was 4.7 and the free-stream pressure about 5 mm 
Hg. The model was at a slight positive angle of attack. In the 
picture the shock-wave pattern (shocks originating from the step 
at the front end of the insert and from the nose cone) is clearly 
and it matches exactly with the shock pattern taken by 
The bright region along 
No measurements 


Fig. 2 
pan film at F 100 sec. 


watts 


visible, 
schlieren pictures without the discharge 
the model seems to be the boundary layer 


FORUM 124 
ee ee ee ee +A i 

ay aj ayn} WN Con AQ 

Baz: 7 Ca 22a 
a," = a,” _ d -— a = f, no” - ie TF Da;;" . é Ea; a 

Fate "+ Ga; F 22b 
Since from Eq. (1), a7 = ar — ag, the terms ay’, ay’, a7’ 
aps”, aps”, age”, and azz” in Eqs. (22a) and (22b) can be rewritten 


as known coefficients using Eqs. (5) through (18 
By equating power coefficients for the wing control surface 
combination the resulting matrix is obtained 
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are available yet to prove this; however, other pictures which 
show typical interaction patterns of shock waves and boundary 
layers tend to confirm it. Regions of dead flow also show a 
brighter discharge. If the current density is increased approxi 
mately five times, striations will appear in some areas of the dis 
charge which seem to indicate the directions of streamlines. As 
the mechanism of the whole phenomenon is rather complicated 
and not well understood, no attempt shall be made here to give 


a physical explanation 
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Hypersonic Flow Characteristics Including 
Real Gas Effects 
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T° A RECENT NOTE! a method has been proposed to extend the 
applicability of the perfect-gas compressible-flow equations 
for air to higher Mach Numbers and stagnation temperatures by 
utilizing suitable approximations for the variation of specific heat 
with temperature. It was also suggested! that the foregoing 
method could be extended for normal and oblique shock proper- 
The purpose of the present note is to point out the existence 
‘ wherein the actual state of the gas, in- 


ties 
of several references? 
cluding real gas effects, has been accounted for in the determi- 
nation of hypersonic equilibrium flow properties and shock char- 
acteristics. 

In reference 2, Feldman, using an IBM 650 computer, presents 
in chart form the equilibrium properties for normal and oblique 
shocks in argon-free air in hypersonic flow and in shock tubes. 
A Mollier diagram based upon National Bureau of Standards data 
and the statistical-mechanical calculations for the speed of sound 
by Logan and Treanor is also included, from which isentropic 
flow properties can be determined. 

In reference 3, a program following the methods of Feldman? 
and Moeckel,® but using the more up-to-date 1956 ARDC Model 


Atmosphere, was run on an IBM 701 computer. The results, 
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hypersonic normal and oblique shock equilibrium properties jg 
cluding the effects of wing sweep, are presented in chart for 
The computations based on the thermodynamic properties caley) 
lated by Hilsenrath and Beckett for argon-free air, which were 
limited to stagnation temperatures above 3,600°R. in referengg 
2, were extended down to 1,440°R. in reference 3 using the calog 
lations in air by Gilmore.’ (Errors on the order of only 1 pe 


cent result from the use of air properties calculated without 
argon.*) 

A similar and somewhat more detailed approach to the caleuk. 
tion of normal shock characteristics in equilibrium air is presented 
by Hochstim in reference 4. Computations were made at re 
duced ambient temperatures and pressures by extrapolating 
NBS atmospheric property data to zero pressure. The results 
which can also be applied directly to oblique shock properties, are 
presented in tabular and chart form 

For the computation of weak oblique shock properties and 
isentropic flow properties, the methods proposed in reference ] 
will furnish a simple and accurate means by which moderate real 
gas effects can be considered and will serve as an effective inter 
mediate between perfect-gas compressible-flow tables and reak 


gas hypersonic-flow charts. 
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